A total of 491 commercial meat rabbits were used in this study, including 205 EWR and 286 NZW. The nutritional levels and feeding management were described in our previous study (Zhang et al., 2011). The growth traits dataset included body weight (BW) of 35 (BW35), 56 (BW56), and 70 (BW70) day of age. Ear tissues were collected and the genomic DNA was extracted using AxyPrep Genomic DNA Miniprep Kit (Axygen, USA), the DNA was stored at −20°C.

To scan the genetic polymorphisms of TBC1D1 gene coding region, five primer pairs were designed according to the TBC1D1 gene mRNA sequence (GenBank accession NO. XM_002709384.1) using Premier 5.0 (Table 1). For each primer pair, 20 rabbits were randomly selected from 491 individuals for mutational analysis. The 30 μL PCR reaction volume included 15 μL 2× Taq PCR MasterMix (TIANGEN, Beijing, China), 3 μL cDNA template (20 ng/μL), 9.6 μL ddH₂O, 1.2 μL of each primer (10 pmol/μL). The PCR reaction was performed with the following conditions: one initial denaturation cycle at 95°C for 4 min, followed by 34 cycles of 95°C for 30 s, 50°C to 61°C for 30 to 45 s and 72°C for 40 to 80 s, end with an extension cycle at 72°C for 10 min. The purified PCR products were directly sequenced on a 3700 DNA sequencer in both directions using the BigDye Terminator sequencing kit (Applied Biosystems, Foster City, CA, USA), according to the manufacturer’s instructions.

Primer pair of HRM214F and HRM214R (Table 1) was designed to amplify the 90 bp fragment containing the target SNP of c.214G>A, which was further subjected to HRM analysis. HRM analysis was performed on the Rotor-Gene Q series real-time PCR detection system (Corbett Research Pty Ltd, Sydney, Australia). The 10 μL reaction volume contained 5 μL 1×SsoFast EvaGreen supermix (Bio-Rad, USA), 1 μL DNA template (20 ng/μL), 3.2 μL ddH₂O, 0.4 μL of each primer (10 pmol/μL), which was pipetted into each well on the 72-well PCR plate. The PCR profile was as following: one initial denaturation cycle at 98°C for 3 min, 35 cycles of 98°C for 5 s, 62.1°C for 30 s, subsequently 95°C for 10 s. Each run contained a non-template control (NTC) and three known genotype individuals including wild-type, heterozygote, homozygote from sequencing verified reference sequence types as the standards. Genotyping of each sample was repeated twice in the same run. The high-resolution melting data was analyzed with the precision melt analysis software (Corbett Research Pty Ltd, Sydney, Australia).

The DNA sequences were assembled and analyzed with the DNAstar program (DNAS Inc, Madison, WI, USA). The localization prediction of rabbit TBC1D1 gene mutations was obtained with SMART web server (Letunic et al., 2009). HRM curve data was analyzed with the Analysis-only Rotor-Gene Q series software VIRTUAL MODE software (Corbett Research Pty Ltd, Sydney, Australia). Allele frequencies, genotype frequencies, heterozygosity (He), effective number of alleles (Ne) and Polymorphic Information Content (PIC) were estimated based on the previous study (Botstein et al., 1980). The x² analysis was used to test for deviation of genotype from Hardy Weinberg Equilibrium (HWE). The associations of genotypes with the growth traits were carried out by the least-squares method as applied in the General Linear Model (GLM) procedure of SAS 9.2 program. The statistical model as follows:

Where Y_ijkl was a record of the trait, μ was the overall mean of observations, G_i was the gender effect, M_j was the fixed genotype effect, P_k was the fixed born time effect, and e_ijkl was the residual error.

Genetic polymorphism of TBC1D1 in rabbit has not been reported before, we first scanned the coding region of the TBC1D1 to detect potential SNPs in rabbit. A total of 12 SNPs were detected and all of them were distributed over five exons region (Figure 1). Among them, there was only one non-synonymous mutation (c.214G>A) which was detected in exon 1 and resulted in the amino acid change from glycine to arginine (G72R). In porcine TBC1D1 gene, three synonymous mutations (g.40A>G, g.151C>T, and g.172G>A) were also observed, which were located in exon 2; other two mutations (g.219G>A and g.252G>A) were located in intron 2 and the g.219G>A mutation was correlated with some growth traits and slaughter traits (Fontanesi et al., 2011). In the human TBC1D1 gene, a non-synonymous mutation (R125W, rs35859249) was observed which showed an association with higher risk of obesity (Stone et al., 2006; Meyre et al., 2008).

The prediction of rabbit TBC1D1 domains showed that the two PTB domains spanned the amino acids 16 to 164 and 168 to 387, respectively; the TBC domain spanned from the amino acids 800 to 1020; a coiled coil spanned amino acids 1077 to 1132, the remains were lower conserved sequences. Among the 12 SNPs, seven synonymous mutations (c.2451G>A, c.2514A>G, c.2535G>A, c.2589G>A, c.2763T>C, c.2844T>C, c.2937G>T) were located in the TBC domain of TBC1D1 in rabbit. The TBC domain is known for its Rab GAP activity and regulation of Rab family small molecular weight GTPases (Will and Gallwitz, 2001). In a previous study, a mutation in the tbc1d1 gene in SJL strain mice resulted in a truncated protein absence of the TBC Rab-GTPase-activative protein domain which suppressed high-fat diet-induced obesity (Chadt et al., 2008). It suggested that the functional mutation in TBC domain may play a potential role in regulating energy metabolism. We detected only seven synonymous mutations in this region; more individuals need to be sequenced to explore the extent of non-synonymous mutations in the TBC domain. The non-synonymous mutation c.214G>A (p.G72R) and other four synonymous mutations (c.153T>C, c.252A>G, c.564G>A, c.591G>A) were located in two PTB domains, respectively. As mentioned previously, the two PTB domains located in N-terminus of TBC1D1 were conserved in AS160. The PTB domain is typically involved in targeting protein to the appropriate cellular location or interacting with a binding partner (Marchler-Bauer et al., 2013). The protein holding PTB domain possesses powerful functions including the ability to bind inositol phosphates and/to other proteins which is dependent on the phosphorylation status of tyrosine residues on the binding partner (Yan et al., 2002). In previous studies, a R125W mutation in TBC1D1 located in PTB domain was associated with obesity susceptibility in humans (Stone et al., 2006; Meyre et al., 2008). Some other studies have proved that altering an R in TBC1D1 at a position corresponding G72 affected the binding properties of the PTB domain (Zhou et al., 1995; Yaich et al., 1998). Thus, we could presume that G72R mutation may be linked to the binding of TBC1D1 to the binding partner known as Akt in a PI3-kinase-dependent manner (Ishiki and Klip, 2005) which would have immediate effect on growth traits. To prove the hypothesis, c.214G>A was selected to carry out the association analysis with growth traits in the rabbit.

The evidence accumulated over the past decade also suggests that synonymous mutations could have functional consequences that can result in an aberrant mRNA splice, affecting mRNA stability and thus protein expression and enzymatic activity. The synonymous mutations were also demonstrated to affect protein conformation and have function and clinical consequence. In future studies regarding TBC1D1 and growth traits, the effects of synonymous mutations should be considered in our research (Sauna and Kimchi-Sarfaty, 2011).

High-resolution melting (HRM) analysis is a relatively new post-PCR analysis method used to identify variations in nucleic acid sequences. This method is based on detecting small differences in PCR melting curves and is enabled by improved dsDNA-binding dyes used in conjunction with real-time PCR instrumentation. In the presence of the fluorescent DNA intercalating dye and unlabeled oligonucleotide probes, the change in stability of DNA duplexes induced by subtle sequence variation is monitored as the temperature of the solution increases (Martino et al., 2010). Therefore, the shape of the melting curve changes with the sequence variation (Montgomery et al., 2007). Since different genotypes have unique melt curves and Tm, it is possible to identify the two homozygotes according to Tm directly and to discern the heterozygotes based on Tm and the changes in the shapes of the melt curves as compared to homozygotes (Nguye-Dumont et al., 2009). In previous studies, genotyping by HRM showed high sensitivity and specificity (Wittwer et al., 2003; Garritano et al., 2009; Martino et al., 2010). To assess the sensitivity and specificity of HRM in the present study, we undertook direct sequencing of 36 samples on the same plate. The sensitivity and specificity of HRM for c.214G>A mutation detection were 0.97 [35 true positive/(35 true positive+1 fasle negative)], and 0.98 [71 true negative/(1 false positive+71 true negative)], respectively. The results were consistent with previous studies. Therefore, HRM analysis was used to genotype the c.214 G>A of TBC1D1 in two rabbit breeds and the normalized graphs produced by analyzing amplified fragment of the three different genotypes are presented in Figure 2.

The genotype and allele frequencies of c.214G>A in exon 1 of TBC1D1 in two rabbit breeds are shown in Table 2. Allele A was the predominant allele with frequencies of 0.78 in EWR and 0.67 in NZW, respectively. The frequencies of genotype AA for EWR and NZW were the highest (0.64 and 0.49, respectively). The polymorphic site showed low heterozygosity (He = 0.35) and moderate polymorphism information content (0.25<PIC<0.50) in EWR, high heterozygosity (He = 0.44) and moderate polymorphism information content (0.25<PIC<0.50) in NZW. The genetic diversity indicated higher genetic variation and a selective potentiality in NZW which could be expected to gain more genetic progress, while the lower genetic variation in EWR was probably due to the long-term artificial selection which made the negative allele frequency reduction (Andersson, 2001; Darwin, 2009). Genotype frequencies of c.214G>A were not in HWE in either EWR (x² = 7.790, p = 0.020) or NZW (x² = 11.002, p = 0.004), the number of AA homozymous subjects was higher than expected in both breeds, the potential cause would be the population dynamics resulted from the inbreeding of livestock that led to an increased frequency of homozygotes (Schaid and Jacobsen, 1999). In order to efficiently analyze genetic variation, more unrelated individuals and breeds need to be included in the study.

In the present study, the association analysis suggested that different genotypes had significant effects on the BW35, rabbits with genotypes GA and AA had higher 35 d BW than those with genotype GG in both EWR (p<0.01) and NEW (p<0.05). But the TBC1D1 genotypes had no significant association with other growth traits (p>0.05) (Table 3). In previous studies, a single nucleotide polymorphism (g.219A) located in intron 2 of porcine TBC1D1 had a significant association with higher ham weight and average daily gain (ADG) in Italian Large White pigs, the less frequent genotype g.219GA showed a lower ham weight and ADG (Fontanesi et al., 2011). In chicken, the TBC1D1 haplotype was significant for body weight at 12 weeks of age in broilers (Rubin et al., 2010; Gu et al., 2011). The TBC1D1 with R125W mutation was associated with obesity susceptibility in childhood and adult humans (Stone et al., 2006; Meyre et al., 2008). These results indicated that the TBC1D1 may be involved in rapid postnatal growth resulted from the energy metabolism in muscle and fat as genotypes at some of SNPs in TBC1D1 were associated with growth traits. Due to the negative genotype GG for c.214 G>A variant in TBC1D1 having a low frequency, groups with genotype GG had lower weight during postnatal growth than groups with AA and GA in both EWR and NZW. These effects indicated that the c.214G>A polymorphism of TBC1D1 was involved in regulating growth traits in rabbit. The groups with allele A had a positive effect on caloric shift which could lead advantageous growth in rabbit.

The proposed function for TBC1D1 was reasonable for its role involved in systemic energy balance by regulating the transport of GLUT4 or GLUT1 into cells. We predicted that the c.214G>A occurring in the protein interaction domain altered interaction of TBC1D1 with a binding partner known as Akt in insulin-induced glucose metabolism. The abnormal interaction may alter caloric disposal in insulin sensitive tissues, furthermore, the genotypes AA and GA had a positive effect on this process. When a high energy source is required for a rapid postnatal growth carried out around weaning, the groups with genotypes AA and GA have higher ability to uptake glucose for predominant binding of TBC1D1 to a binding partner and accumulate higher WB35 in two breeds. This speculation is supported in New Zealand obese (NZO) mice which present a feature of early-onset obesity and are highly susceptible to weight gain when fed a high-fat diet. Studies regarding NZO identified some dietary fat-responsive genetic loci and TBC1D1 was one of the genes in the Nob1.24 segment which was a suppressor of high-fat diet-induced obesity and the TBC1D1 gene mRNA expression was upregulated in high-fat diet energy metabolism (Crofford and Davis, 1965; West et al., 1992; Wuschke et al., 2006). Therefore, when a well-distributed growth is carried out in early development before 70 d of age for rabbit, the binding of TBC1D1 to Rab targets in glucose metabolism induced by energy source keeps a relatively stable state that had no significant effect on WB56 and WB70 in the present study. These results were also consistent with the function of TBC1D1 that is a regulator of either glucose or fatty acid transporter to achieve systemic energy balance by regulating the calories in cells (Kahn and Klip, 2000; Raz et al., 2005). Obviously, further experiments would be carried out to prove the precise molecular mechanism of TBC1D1 as well as the PTB domains of TBC1D1 involving in insulin-stimulated glucose metabolism, including the comparison of enzymatic activity among the alternative genotypes.