1 INTRODUCTION

Meretrix petechialis is one of the most important commercial clam species farmed in China (Liu et al., 2006). Despite continuous progress and improvement in artificial reproductive technologies and aquaculture techniques, clam diseases remain a major limiting factor in the clam culture industry. Outbreaks of vibriosis in clams often affect the profitability and sustainability of aquaculture operations. From a sustainability perspective, selective breeding for disease resistant clams is an attractive strategy for disease prevention (Liang et al., 2017). However, in contrast to many other economically important traits (e.g., growth rate, feed conversion efficiency), genetic improvement of disease resistance has been hindered by the difficulty of measuring appropriate phenotypes accurately (Bishop and Woolliams, 2014). Therefore, markers accurately reflecting disease resistance and having a significant genetic component are necessary in resistance selection.

To date, some association analyses between molecular markers and Vibrio resistance traits have been reported in M. petechialis (Nie et al., 2013, 2015; Zou and Liu, 2016). However, the effectiveness of these markers is minor, accounting for only a small fraction of the genetic variation (Yang et al., 2010). Thus, more attempts are needed to develop markers that are closely linked to Vibrio resistance. In previous studies, the use of gene expression profiles as indirect tests for improved disease resistance has been proposed for selective breeding in Atlantic salmon (Robinson and Hayes, 2008). Research on M. petechialis has demonstrated that hundreds of genes are involved in the physiological and immunological response to bacterial infection (Jiang et al., 2017). Although numerous studies have demonstrated differences in gene transcription among genetic groups or lines in aquatic species (Langevin et al., 2012; Marancik et al., 2014; Reyes-López et al., 2015; Robledo et al., 2016), few have partitioned genetic variation and estimated heritability for gene transcription (Roberge et al., 2007; Normandeau et al., 2009; Aykanat et al., 2012; He et al., 2017).

Gene expression variation can be considered as a phenotype (Houle et al., 2010), and as with any traditional phenotypic trait, gene expression variation is influenced by various genetic and environmental factors (Leder et al., 2015). It is clear that gene expression variation is widespread among individuals and taxa and can be heritable (Whitehead and Crawford, 2006; Dixon et al., 2007; Powell et al., 2013). Additive genetic effects of gene expression represent evolutionary potential, which allows us to predict the magnitude of its response to selective pressure (Leder et al., 2015; He et al., 2017). Nonadditive effects also contribute significantly to phenotypic variation but are not heritable (Evans and Neff, 2009). All these components of genetic parameters are important contributors to phenotypic variation. Therefore, it is fundamentally important to disentangle the contributions of additive genetic and non-additive genetic variance to overall gene expression variance. This information of genetic parameter is critical for evaluating the potential for gene-expression performance markers to be used in selection programs targeting resistant clam broodstocks.

Recent studies have shown that the basal expression levels of one set of Vibrio resistance-related genes (Big-Def, BIRC7, NFIL3, CTL9, and Bax) significantly correlated with survival in challenge tests in clam M. petechialis (Jiang et al., 2018). Therefore, it seemed that selection decisions could be made based on the expression level of resistance-related genes that are genetically correlated with disease resistance. The indirect selection by gene expression trait would simplify the data collection, avoid the challenge testing, and reducing costs and time consumption. Theoretically, higher additive genetic effects reflect higher adaptive potential in response to selection pressures; thus, the gene transcription traits closely linked to phenotype and with higher additive genetic component are better candidates for an indirect selection criterion in practice. However, the genetic variance components for the transcription of resistancerelated genes and their genetic correlations with body size have not been evaluated in M. petechialis. Therefore, in this study, an animal model was used to estimate the genetic variation component for Vibrio resistance-related genes, and estimate the genetic correlations between Vibrio resistance-related genes and body sizes. The results could provide some key parameters for understanding of the potential resistant selection response of gene transcription for resistancerelated genes in non-model organisms.

2 MATERIAL AND METHOD 2.1 Experimental families

Experimental families were produced in a hatchery at the Zhejiang Mariculture Research Institute (Wenzhou, China). In July 2014, 11 two-year-old males and 11 two-year-old females were randomly selected from the clam population and used to produce 23 full-sibling families according to a slightly modified Berg and Henryon design (Berg and Henryon, 1999). For example, each sire was mated with 2–3 dams, and each dam was mated with 2–3 sires. Thus 11 dam halfsib families and 11 sire half-sib families were produced, the detailed breeding design for these fullsibling families was presented in Fig. 1. Detailed rearing conditions for larvae have been described by Wang et al. (2011). The full-sibling families of juvenile clams were randomly distributed in pond at the Zhejiang Mariculture Research Institute. In the farming process, different families were reared separately by mesh in the same pond and exposed to similar environmental conditions to minimize the effects of uncontrolled environmental factors. In detail, the stocking density of each family was the same, and the average survival rate of families was up to 92.4%±1.9% at 15-month-old, which means that no accidental mass mortalities occurred in the cultivation.

2.2 Sample collection

Before the estimation of genetic parameters, six tissues (digestive gland, mantle, gill, foot, adductor muscle, and hemocytes) were collected from four healthy clams for analysis of the expression of the detected genes between tissues. Hemolymph (~500 μL) was withdrawn from each adult clam with a 1-mL disposable syringe, and was then centrifuged at 3 000×g for 10 min at 4℃. The pellet was used in total RNA extraction. Another five tissues from each clam were taken, dissected, immediately frozen in liquid nitrogen, and stored at -80℃ until tissue distribution analysis could be conducted.

The 12-month- and 15-month-old clams from two batches were used to estimate genetic parameters respectively. We randomly selected five individuals from each family from each sampling batch for the genetic analysis. The shell width (SW) and total weight (TW) of 12-month-old clams were measured using a digital Vernier caliper and a digital balance respectively. After measuring the SW and TW, the digestive gland from each clam was dissected, immediately frozen in liquid nitrogen, and stored at -80℃ until further processing.

2.3 RNA extraction and cDNA synthesis

Total RNA was extracted from samples using an EZNA^® Total RNA Kit II (Omega Bio-Teck, USA) according to the manufacturer's instructions. RNA degradation and contamination were monitored on 1% agarose gels and then quantified using a NanoDrop-1000 spectrophotometer (Thermo Scientific, USA). cDNA was synthesized from total RNA (1 μg) using a PrimeScriptTM RT reagent kit with gDNA Eraser (TaKaRa, Japan) according to the manufacturer's protocol.

2.4 Semi-quantitative RT-PCR

Semi-quantitative RT-PCR was performed to determine the expression levels of Vibrio resistancerelated genes in different tissues from adult clams. The elongation factor 1α gene (EF1α) was used as the internal reference gene in semi-quantitative RT-PCR analysis (Jiang et al., 2017). Primer sequences are listed in Table 1. PCRs were performed in a 25 μL reaction volume containing 1 μL cDNA, 0.2 mmol/L of each primer, and 12.5 μL Premix Taq^TM (TaKaRa, Japan). The appropriate number of cycles of these six genes were selected, 18 cycles for EF1α, 26 cycles for Bax, NFIL3 and Big-Def and 28 cycles for BIRC7 and CTL9. These cycling parameters could ensure that the amplification product is clearly visible and can be quantified. That is the amplification is in the exponential range and has not reached a plateau yet. The cycling protocol was one cycle of 94℃ for 4 min; 18, 26 or 28 cycles of 94℃ for 15 s; 55℃ for 15 s; and 72℃ for 15 s, and the final extension time was 7 min. PCR was conducted using a Takara thermocycler. PCR products were separated via electrophoresis in 2% agarose gels that were stained with 4S Red Plus Nucleic Acid Stain (Sangon, China) and visualized under ultraviolet light. The quantification of bands was performed by Bio-Rad Quantity One.

2.5 Quantitative real-time PCR

Based on a previous study (Jiang et al., 2018), 5 Vibrio resistance-related genes were selected for detection of the genetic parameters of gene expression levels using clams randomly selected from the 23 full sibling families by quantitative real-time PCR (qRTPCR). The β-actin and EF1α genes were used as internal references to normalize the relative expression levels between samples (Jiang et al., 2017). The primers for qRT-PCR are the same used for semiquantitative RT-PCR. qRT-PCR was performed in QuantStudio 6 Flex (Applied Biosystems, USA) machine using a QuantiNova SYBR Green PCR Kit (Qiagen, Germany). Each reaction was performed in triplicate under the following conditions: 95℃ for 2 min, 40 cycles of 95℃ for 5 s, and 60℃ for 20 s, followed by melting curve determination. The ΔCt value was used to quantify the relative gene expression value, the greater the ΔCt value, the lower the gene relative expression level (Pfaffl, 2001).

2.6 Estimation of genetic parameters

The phenotypic variance of 12-month- and 15-month-old clams were estimated using the following animal linear model:

where Y_ijk is the phenotypic value of k^th offspring of i^th sire and j^th dam; μ is the mean phenotype of the sample; a_k is random additive genetic effect of clam k; I_ij is the no additive effect due to the interaction of i^th sire and j^th dam; and e_ijk is the residual of the k^th individual of i^th sire and j^th dam.

For the animal linear model, heritability was calculated as follows:

where σ_a² is additive genetic variance components; σ_I² interaction between sire and dam, also named family variance components; and σ_e² residual variance components.

Gene transcription variance was partitioned into additive genetic variance (σ_a²) and the interaction between sire and dam (σ_I²) separately. Wilcoxon signed-rank test was conducted to compare the genetic parameters of gene transcription trait between 12-month- and 15-month-old clams.

2.7 Estimation of phenotypic and genetic correlations

The phenotypic variance SW and TW were partitioned using the following animal linear model: and the parameters were as described above. A bivariate linear animal model was used to estimate the variance and co-variance components for the five transcription traits and SW and TW. The phenotypic or genetic correlation between each paired trait was calculated as: where σ₁₂ is the phenotypic or genetic covariance between two traits, σ₁² and σ₂² are the additive phenotypic or genetic variances of trait 1 and trait 2, respectively. All the data were analyzed using the ASReml 3.0 software package in R (Gilmour et al., 2009).

3 RESULT 3.1 Tissue expression profile analysis

The expression levels of the five genes varied among different tissues (Fig. 2). For example, CTL9 had the highest expression level in the digestive gland, relatively low expressive abundance in hemocytes and adductor muscles, and no expression in the mantle, gill, and foot. The other four genes were expressed in all six tissues, but expression levels varied among tissues. For example, the apoptosisrelated genes Bax and BIRC7 both had the highest expression levels in the gills but were only moderately high in the digestive gland. The transcription factor of inflammatory response NFIL3 had high expression level in the mantle, gill, and foot and relatively low expression in the other three tissues. Big-Def had low expression in the gills and high expression levels in other five tissues, with the highest in the digestive gland. These results revealed that expression of all five genes could be detected in the digestive gland and hemocytes in M. petechialis. Considering that few hemocytes could be extracted from the clams, the digestive gland was selected as the major source for the following gene expression analysis.

3.2 Characteristic of the gene transcription traits

The number of observations, means, standard deviations, coefficients of variation (CV) for the five gene transcription traits (Bax, BIRC7, Big-Def, CTL9, and NFIL3) in two growth stages are shown in Table 2. The expression level differed among the five genes in digestive gland. Big-Def had the highest relative expression level and the highest CV in 12- and 15-month-old clam group, whereas the relative expression level of BIRC7 was the lowest. Furthermore, the gene expression levels of 12- and 15-month-old clam were compared. With the exception of BIRC7, other genes showed lower expression in 15-month-old clam group. The gene expression difference for five genes between 12- and 15-month-old clam group was 8%–23%. Compared to other genes, the expression of BIRC7 and NFIL3 for 12-month-old was more similar to that of 15-month-old clam group. The result of t-test suggests that expression level of the 5 genes between 12- and 15-month old clam were statistically and significantly different (P < 0.05).

3.3 Genetic estimates of the transcription traits

Estimation of the variance components and heritability of gene transcription trait in 12-monthand 15-month-old clams are presented in Table 3. The gene transcription trait of two different growth stages exhibited substantial similarities in genetic variance components. Across all five genes, the gene transcription traits of two stages were no significant different in additive genetic effects (Wilcoxon signedrank test, P=0.31) and sire-by-dam interaction effects (no additive genetic effects) (P=0.13). In detail, the BIRC7 gene showed significant additive genetic variations in the 12-month and 15-month-old clams, correspondingly the high heritabilities 0.84±0.32 and 0.91±0.34, respectively. Whereas, the heritabilities of other four genes (Big-Def, BAX, NFIL3, and CTL9) were low-to-moderate but not significantly. The NFIL3 gene showed significant sire-by-dam interaction effects.

3.4 Genetic correlation for transcription trait and body size

The phenotypic and genetic correlations were analyzed based on the body size SW, TW and expression of BIRC7 gene in the stage of 12-monthold clams. The phenotypic and genetic correlation coefficients between SW and BIRC7 gene were 0.05±0.15, -0.04±0.44, the phenotypic and genetic correlation coefficients between TW and BIRC7 gene were -0.18±0.15 and -0.97±0.98 respectively. No significant phenotypic and genetic correlations were detected between BIRC7 gene transcription trait and body size in the 23 full sibling families.

4 DISCUSSION

The basal gene expression level could vary markedly in different organs of the same animal (Engwerda and Kaye, 2000). Therefore, the choice of research organ is crucial for gene expression studies. In our experiment, all the five resistance-related genes were highly expressed in hemocytes and digestive gland in M. petechialis. It has been shown that hemocytes play a key role in animal immune response (Lavine and Strand, 2002). However, in M. petechialis, due to few hemocytes extracted from the clams' bodies, it is not a preferred organ at the individual level or for batch analyses. In mollusks, the digestive gland is another important tissue of the immune system because their epithelial cells can produce major immune molecules (Röszer, 2014). In addition, the digestive gland has been used as a preferred research organ in M. petechialis (Gao et al., 2012; Zou and Liu, 2016). In this study, the digestive gland is the preferred tissue for gene expression analysis considering that all five genes expressed in digestive gland and it could be collected conveniently.

Gene expression variation can underlie complex phenotype variations (Ayroles et al., 2009) and play important roles in the response to environmental change and stress (Crawford and Powers, 1992). It has long been thought that a large part of gene expression variation is highly heritable with an underlying polygenic architecture (Cheung and Spielman, 2002; Gilad et al., 2008). In recent years, the genetic parameter of gene transcription has been estimated in some aquatic species (Aykanat et al., 2012; Brokordt et al., 2015; Leder et al., 2015; He et al., 2017). Additionally, the heritability of different gene transcription was diverse. A potential explanation is that different phenotypic or transcription traits undergo different strength of selection throughout an animal's life. For example, as those closely related with the biological fitness, i.e. for survival and reproduction etc., undergo more strong selective pressures than other traits, and usually show lower heritabilities. Herein, we examined the genetic parameters of gene transcription for five Vibrio resistance-related genes in full sibling clam families. The results showed that the heritability of BIRC7 was 0.91 and the heritability of other four genes was not different from zero, which showed that different gene transcription traits could exhibit dramatic variation in heritability. In previous studies, Aykanat et al. (2012) found that the heritability of four cytokine genes in Chinook salmon varied from 0.05 to 0.44, and Tedeschi et al. (Tedeschi et al., 2016) reported a higher heritability of three heat shock protein genes in sea turtles (mean h²=0.58). This presence of heritability suggests that Vibrio resistance-related gene expression variation of BIRC7 gene are heritable. Thus, BIRC7 gene transcription may provide reliable genetically based markers for use in selective breeding and the performance of gene transcription can be improved by successive selection.

Genetic analysis is heavily dependent on data, a larger data set will provide a more accurate estimation of genetic parameters, but the costs used in getting the data will also be increased. No research yet has provided the evidence that appropriate sample size is required for genetic analysis. Due to the measuring procedure is complicated for transcription traits and the high cost, the sample size used in genetic analysis is usually not large. For example, Aykanat et al. (2012) used a data set including 48 families and 192 individuals to partition the transcription variation of four cytokine genes into additive genetic, nonadditive genetic, and maternal components in Chinook salmon Oncorhynchus tshawytscha. In our study, we used 23 full sibling families to estimate the genetic parameters of 12-month- and 15-month-old clams' transcription traits. The results showed that there was no significant differences in additive genetic effects and sire-by-dam interaction effects between two growth stages (Wilcoxon signed-rank test P=0.31, P=0.13), which indicated that the genetic variance components of gene expression traits were relatively stable during short-term growth. However, the estimation of heritability in our study was accompanied by high standard error, which are likely associated with low sample sizes, and possibly resulted in our failure to detect significant low-to-moderate heritabilities (Table 3). Nevertheless, we observed significant additive effects in BIRC7, which might be valuable in the genetic improvement of Vibrio resistance in artificial selection approaches.

Animals are frequently faced with trade-offs between growth and mortality rates. Enhancement of disease resistance may become with compromise in growth traits. For instance, Yáñez et al. (2016) has reported a significant negative genetic correlation between resistance against Piscirickettsia salmonis and harvest weight in coho salmon. Furthermore, a slightly negative genetic correlation between resistance to vibriosis and harvest weight was found in Atlantic cod (Bangera et al., 2011). These results suggest that selective breeding for faster growth may have a negative effect on disease resistance in these species. Conversely, a significant positive genetic correlation between resistance to fish pasteurellosis and body length was reported in the gilthead sea bream (Antonello et al., 2009), which may be due to resistant individuals experiencing superior growth. Additionally, genetic correlations that are not significantly different from zero have also been found between disease resistance and harvest weight in some aquatic species (Silverstein et al., 2009; FloresMara et al., 2017). Results of these studies suggest that genetic correlations between disease resistance and growth of different species are not uniform. Therefore, genetic correlation must be evaluated case by case in practical breeding. In this study, we found that there is no significant genetic correlation between Vibrio resistance-related gene BIRC7 and body sizes. Our results indicate that variation in gene transcription could help the genetic improvement of Vibrio resistance and will not influence body size in M. petechialis.

5 CONCLUSION

In conclusion, our results first demonstrate the variation in the expression of the Vibrio resistancerelated BIRC7 gene in M. petechialis is heritable, indicating the feasibility of improving this trait by artificial selection. Moreover, no significant genetic correlation is detected between the Vibrio resistancerelated BIRC7 gene and body sizes, suggesting that artificial selection for Vibrio-resistance will not influence body size in M. petechialis.

6 DATA AVAILABILITY STATEMENT

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study. All data generated and/or analyzed during this study are included in this published article.