Appropriate growth models to describe early growth of Kejobong goat based on Growth Hormone ( GH ) gene sequence analysis

The objectives of this study were to reveal appropriate growth models describing early growth of Kejobong goat based on Growth Hormone (GH) gene sequence analysis. A total of 35 DNA samples and 1.960 records of quantitative traits of Kejobong goat were collected. The exon 3 of GH gene was amplified and was sequenced to determine the SNP. Body weight and body measurements of the goats were taken at 0-14 weeks of age. Four non-linear growth models were applied for analysis of growth to compare growth performance of different genotypes by Non-Linear Mixed Model. A non-synonymous mutation (g1170AG) genotyped into GG, AG and AA was significantly associated with growth traits. Animals with heterozygous genotype AG showed higher growth traits than animals with homozygous genotype AA. Nonetheless, animals with homozygous genotype GG had the same growth traits with those animals with heterozygous genotype AG and homozygous genotype AA. The most fitted model for describing body weight was Von Bertalanffy model, while for describing wither height and hip height GH gene and growth model analysis on goat (Sutopo et al.) 124 was Brody model. SNP at exon 3 of the GH gene can be used as genetic marker for improvement of growth traits of Kejobong goats.


INTRODUCTION
Kejobong goat is one of indigenous Indonesian breeds, which only exists in Purbalingga District, Central Java Province, Indonesia, and it is conventionally raised by local farmers. This goat belongs to Southeast Asian lineage and is confirmed to be descendant from crossbred of Kacang and Etawah Grade goats (Kurnianto et al., 2012;Kurnianto et al., 2013;Lestari et al., 2018 B ). As the meat animals, Kejobong goat had 41.30% carcass yield that comprised 67.06% meat and 32.94% bone, while its meat is known to have less cholesterol than meat of Kacang and Etawah Grade goat (Aqsa et al., 2011). Kejobong goat is popular at the district because of its high rate of growth, good reproductive performance, high resistance to local diseases and parasites and ability to survive and growing ability under poor feeding conditions (Kurnianto et al., 2012;Febriana et al., 2017). However, the breeders often have difficulty to satisfy the market demand on slaughtering weight. This is probably due to limited information of appropriate breeding strategy for Kejobong goat to accomplish the breeding goal of high meat productivity.
Growth analysis can provide valuable information about mature weight, growth rate and mature time. Growth rate and body weight of animal at different ages influence productivity of meat and have deterministic effects on the profitability of meat production (Kheirabadi and Rashidi, 2019). Particularly, growth rate has large effect on meat producing efficiency up to slaughtering age which is crucial for economic success of animal production (Abbasi et al., 2012). According to Junior et al. (2013) and Ripoll et al. (2016), animals that have a large frame size of body tend to have higher potential of growth and have a higher proportion of meat. Therefore, besides body weight, body size is also important trait to be considered for performing animal selection. Study of growth analysis has been done by previous researchers (Waheed et al., 2011;Setiaji et al., 2013;Zadeh et al., 2015;Raji et al., 2015;Lupi et al., 2016;Zadeh and Gorbani, 2018;Ghiasi et al., 2018;Rout et al., 2018;Kheirabadi and Rashidi, 2019), however they were only using phenotypic data into analysis. In this study, conventional growth analysis was modified by including genotype records to growth analysis.
Early growth of kids is an economically important trait that affecting profitability in goat production (Baranzadeh et al., 2012;Moghbeli et al., 2013;Sadeghi et al., 2019). Physiologically, growth of an animal is a result from a complex process of metabolism including a coordinated action of several hormones that controlled by expression of their responsible genes (Mahrous et al., 2018). Growth Hormone (GH) gene is one of numerous genes which have large effect on growth performance of an animal. GH gene is encoding growth hormone that produces in anterior pituitary and is necessary for postnatal growth and metabolism in vertebrates (Ge et al., 2003). This hormone is known to have a broad impact on biological activity in all body cells, such as controlling and coordinating the flow rate of metabolic process, enhancing glycogen, protein, DNA and RNA biosynthesis and promoting the deposition of fat and the disintegration of fatty acids and glucose in the tissue (Gorlov et al., 2017;Wickramaratne et al., 2010;Othman et al., 2015;Seevagan et al., 2015;Singh et al., 2015). Therefore, GH gene is considered to be a prime factor which affects growth performance of an animal.
Based on these backgrounds, effect of GH gene on growth traits, especially from a point of genetic improvement is important to build breeding plan for high meat productivity. Prospectively, result of this study is not only suggesting appropriate management practice for improving production for the breeders, but also providing information of genetic marker in Kejobong goat for breeding selection in the future through Marker-Assisted Selection (MAS) and/or Marker-Assisted Introgression (MAI) and appropriate mathematical growth models of Kejobong goat. Therefore, the objective of this study was to reveal appropriate growth models describing early growth of body weight and body measurement of Kejobong goat based on the effect of Growth Hormone (GH) gene sequence analysis.

Ethical approval
All procedure involving animals were based on the standard rule of animal treating as appointed in the Republic of Indonesia's law, number 41, 2014.

DNA extraction, Polymorphism Chain Reaction (PCR) and sequencing
Blood samples for DNA analysis were taken by 3cc spuit from jugular venous that previously cleaned with alcohol. The blood was then collected in vacutainer blood collection tubes with an anticoagulant (EDTA). DNA was extracted from whole blood by gSYNC DNA mini kit (Geneaid Biotech, Taiwan) according to the manufacturer's standard protocol for PCR and sequencing analysis. Forward primer F: 5'-TAGAAATGGGGGTGTGTGGGGT-3' and reverse primer R: 5'-CATCCTCCACTGCCATCCAACA-3' (Sigma-Aldrich, Japan) were used to amplify GH gene exon 3. PCR was carried out in total volume 50 μL comprising 1 μL KOD Plus (Toyobo, Japan), 5 μL buffer, 5 μL dNTP, 2 μL MgSO4, 1.5 μL forward primer, 1.5 μL reverse primer, 32 μL PCR water and 2 μL DNA template. PCR amplification was running with an initial denaturation at 94°C for 2 min, followed by 40 cycles of denaturation at 94°C for 15 sec, primers annealing 66.8°C for 30 sec and extension at 68°C for 19 sec. PCR products were electrophoresed using 1.3% Agarose gel at 110 V for 20 min. PCR products were then visualized by UV trans-illuminator and was sequenced through Fasmac sequencing service, Japan.

Data analysis
Allelic and genotypic frequencies were directly calculated. Hardy-Weinberg Equilibrium (HWE) was tested using chi-square statistic (χ 2 ) as follows: where χ 2 is the Chi square value; oi the observed value of genotype frequency, ei the expected value of genotype frequency, χ 2 the table using 5% significance level for HWE test.
Heterozygosity (H) was calculated as follows: where H is the value of heterozygosity and pi the frequency of the i th of k alleles.
Sequencing result alignment was analyzed by Clustal W (Thompson et al., 1994) with Molecular Evolutionary Genetics Analysis (MEGA6.0) (Tamura et al., 2013) to find out the SNP within animals. Sequencing result then was translated into amino acids form by standard genetic code to identify amino acid alteration that caused by SNP.
Linear Mixed Model (LMM) was used to analyze association between genotype with quantitative traits by MIXED procedure in Statistical Analysis System (SAS 9.3) (SAS Institute Inc, 2011). The model was: yijkl = µ + Gi + Fj + uk + b1ɑijkl + b2ɑ 2 ijkl + eijkl, where yijkl is the observed value of a dependent variable (body weight or body measurements); µ the overall mean of the population; Gi the fixed effect of i th genotype (i = 1 for GG, 2 for AG, 3 for AA); Fj the fixed effect of j th farm group (j = 1, 2, 3, 4); uk the random effect of k th animal; b1 and b2 the linear and quadratic coefficients of partial regression, respectively; l th individual measurement, ɑijkl age in days of a covariate and eijkl the random residual for yijkl. Difference in the least square means of the genotypes was tested by the Tukey-Kramer (Tukey, 1949).
The nonlinear growth models comprised Brody (Brody, 1945), Von Bertalanffy (Bertalanffy, 1938), Logistic (Verhulst, 1838) and Gompertz (Gompertz, 1825) and they were compared by describing animal growth (Table 1). Growth models were analyzed using Nonlinear Mixed Model (NLMM) by NLMIXED procedure of SAS 9.3 (SAS Institute Inc, 2011). Body weight or body measurements as dependent variables are influenced by genotype and age. Therefore, dummy variables were created to assess the effect of qualitative variables on dependent variables according to the method by Filho et al. (2014). The NLMIXED procedure was used in this study due to its flexibility in engaging the variance covariance structure which could not be identified by traditional regression approach. Also NLMIXED has ability to handle unbalance data (Aggrey, 2009;Galeano-Vasco et al., 2014). This procedure can reduce potential biases despite selective sampling and supply supplemental parameters that characterize variation between individual animals (Craig and Schinkel, 2001).
The models were tested for goodness of fit using -2 log likelihood, Akaike Information Criterion (AIC) (Akaike, 1974), Bayesian Information Criterion (BIC) (Schwarz, 1978) and the residual variances (σ 2 e). AIC and BIC were calculated by the following formula: where n is the number of observation; SSE the Sum Square Errors and k the number of parameters. Smaller values of AIC, BIC or σ 2 e indicate the best fit of the model to the observations.

RESULTS
Result showed that a total 117 bp of GH gene exon 3 encoding 38 amino acid sequence were well amplified. Sequencing result revealed 5 SNPs as transition mutation in parsimonious form, which were g1121AG (SNP1), g1148TC (SNP2), g1160AG (SNP3), g1170AG (SNP4) and g1178CT (SNP5). Genotype frequencies of Kejobong goats were not different from HWE, and the frequency of heterozygosity was 49% ( Table 2). The estimated allele of the GH gene exon 3 in this study was 57% and 43% for G and A, respectively.
Frequencies of genotypes GG, AG and GC were 37%, 40% and 23%, respectively, so that G allele and heterozygous genotype AG were predominant in this locus.
Test of significance showed that the fixed effect of genotype together with effect of farm and linear and quadratic coefficients of age were statistically significant (P<0.05) in BW while the fixed effect of genotype and linear and quadratic coefficients of age were statistically significant (P<0.05) in WH, and the fixed effect of genotype, age of doe and linear and quadratic coefficients of age were statistically significant (P<0.05) in HH. Conversely, the fixed effect of genotype was not significant in CD, CW, HW and HG (Table 3).
Animals of genotype AG demonstrated the highest BW, HW and HH then it followed with animals of genotype GG and AA (Table 4). Comparing genotypes at different periods, BW0 and BW6 in animals of genotype AG (4.35 kg and 8.40 kg) were significantly heavier (P<0.05) than  animals of genotype AA (3.40 kg and 6.59 kg), while animals of genotype GG (3.84 kg and 6.97 kg) showed no significant difference with genotype AG and AA. However, there were no significant effect of genotype at BW2, BW4, BW8, BW10, BW12 and BW14. Significant difference between genotypes for body measurements were observed in wither height (WH6, WH10, and WH14) and hip height (HH12 and HH14). Similar to body weight, animals of genotype AG had significantly (P<0.05) higher wither and hip heights than those animals of genotype AA, but there was no significant difference between animals of genotype GG with animals of genotype AG and AA (Table 5).
Estimated parameters for body weight, wither height and hip height are presented in Table 6, respectively. Growth analysis showed that Von Bertalanffy model had the lowest -2 log likelihood, and two criteria AIC and BIC compared with the other models indicating this model as the best model for describing growth of body weight in Kejobong goat. On the other hand, GH gene and growth model analysis on goat (Sutopo et al.) 128 the highest -2 log likelihood, AIC and BIC were obtained in Logistic model. Brody model in this study showed fit to wither height and hip height well according to its value of -2 log likelihood, AIC and BIC, which was lower than Gompertz, Logistic and Von Bertalanffy model. Von Bertalanffy model fitted best to body weight, estimated 23.01 kg mature body weight (a), 0.39 integration constant (b) and 0.01389 growth rate (k). The best estimated for wither height and hip height by Brody model were 54.65 cm and 58.91 cm for parameter a; 0.37 and 0.36 for parameter b; 0.01577 and 0.01647 for parameter k. In this study, the estimated parameter b for body weight, wither height and hip height were 0.39, 0.37 and 0.36 respectively.
Furthermore, estimated parameter k of body weight in Von Bertalanffy model was 0.01389. Negative correlation was found between parameter k and parameter a (Table 7). This result was confirmed by the fact that Brody model in this study had the slowest parameter k (0.006948) in body weight, yet it had the highest estimated parameter a (25.29 kg) among the others. Similarly, the highest parameter k in wither height (0.2606) and hip height (0.02341) had the lowest estimated parameter a (50.78 cm, 56.26 cm) in Von Bertalanffy and Logistic models respectively.
Representing variability among individual animals, estimated animal variance (σ 2 u) of the body weight in this study was 5.78. The higher the variance, the greater the difference is realized among animals. Furthermore, residual variance (σ 2 e) of the body weight in this study was 0.31 that indicated the gap between predicted value and observed value. Repeatability of body weight by intra-class correlation was 0.95 this study.

DISCUSSION
SNP2 of this study was also found in a report analyzing GH gene of Chinese goat . Among five SNPs in this study, translation result showed SNP4 causing amino acid alteration which changes amino acid sequence in GH gene exon 3. SNP4 as non-synonymous mutation changed the first triplet codon of AGC encoding Serin to GGC encoding Glycine (Figure 1) and we used it to distinguished as GG, AG and AA genotypes, whereas the other SNPs were silent mutation (SNP1 CAG>CAA(Gln); SNP2 TCT>TCC(Ser); SNP3 CCA>CCG(Pro); SNP5 AAC>AAT(Asn)). According to Nei and Kumar (2000), most of synonymous amino acid was found due to substitution of nucleotides in the third codon, while substitution of nucleotides in the first and second codon generate non-synonymous amino acid. Therefore, nonsynonimous mutation that change amino acid sequence in exon region may change the peptide sequence of the encoded protein and influence the function of the protein, which was growth hormone in this study. This hormone has substantial metabolic effects on somatic growth, stimulation of protein synthesis and cellular uptake of amino acids and development of body composition (Hjortebjerg et al., 2017).
Result in this study agreed with a result by Dayal et al. (2016), in which goats with heterozygous genotype AC had the heaviest body weight among five observed genotypes in Black Bengal goat. Gorlov et al. (2017) reported their study in Salsk sheep that sheep with AB genotype significantly had heavier body weight, average daily gain and carcass weight than sheep with AA genotype. A contradictory result was GH gene and growth model analysis on goat (Sutopo et al.) 130  reported by An et al. (2011) that goats with homozygous genotype AA significantly had higher body weight than those of heterozygous genotype AB at age of one and three months old in Chinese goat, however, wither height showed no significant difference. The different results seem to be due to genetic difference that leads to different structure of GH gene and limited number of observations. Therefore, further study is necessary to validate the predominant effect of heterozygote of GH gene with a larger number of animals and more sampled observations. The best model for describing growth of body weight in this study was different with previous study by Kheirabadi and Rashidi (2019), reported that Logistic model fitted worst to body weight, while Brody model fitted most accurately to body weight of Markhoz goat. In this study, estimated mature body weight (a) was 23.01 kg implying that Kejobong goat had heavier mature body weight than Raeini Cashmere goat (17.97 kg) (Ghiasi et al., 2018) and Nondescipt goat (6.42 to 10.55 kg) (Raji et al., 2015) but lighter body weight than Beetal goat (23.39 kg) (Waheed et al., 2011). Those values of parameter b for body weight, wither height and hip height in this study were described to represent the proportion of mature weight attained after birth, calculated by the initial weight and age value (Lupi et al., 2016). On the other hand, Ghiasi et al. (2018) stated that parameter b is a scale parameter that has no biological interpretation. Waheed et al. (2011)   2018) and 0.0108 in Repartida goat (Pires et al., 2017) by applying Gompertz model, so that Kejobong goat in this study is considered to attain mature weight later than Beetal and Cashmere goats but earlier than Repartida goat. Negative correlation between parameter k and parameter a indicated the slower growth rate, the larger mature weight, vice versa. Previous studies supported this results. Kurnianto et al. (1998) reported that animals with slower growth rate tended to have estimated heavy body weight at maturity. Brown et al. (1976) stated that selection for increasing growth rate tended to decrease mature weight, yet its antagonistic association could be minimized by cross-breeding and improving feed quality. On other hand, previous study by Ghiasi et al. (2018) showed lower animal variance (1.29) and higher residual variance (8.01) on growth analysis of Raeni Cashmere goat using Gompertz model than the present study. Repeatability of body weight in this study was higher than repeatability value of South African Angora goat (Snyman and Olivier, 1999) and Boerawa goat (Beyleto et al., 2010). The high repeatability in this study may be resulted by the fact that systematic factors affecting body weight were fitted as many as possible in NLMM and that earlier body weight was a component of latter body weight.
The NLMIXED procedure used in this study has flexibility in engaging the variance covariance structure which could not be identified by traditional regression approach. Also NLMIXED has ability to handle unbalance data (Aggrey, 2009;Galeano-Vasco et al., 2014). This procedure can reduce potential biases despite selective sampling and supply supplemental parameters that characterize variation between individual animals (Craig and Schinkel, 2001). Therefore, this procedure can facilitate growth analysis by including genotype information and estimates accurately the growth performance of Kejobong goats.

CONCLUSION
SNP g1170AG in GH gene is associated with growth traits and can be used as genetic marker for animal selection to improve goat's growth performance. Animals with heterozygous genotype AG showed higher growth performance than homozygous genotype AA. Nonetheless, animals with homozygous genotype GG showed no difference with either heterozygous genotype AG or homozygous genotype AA. Model (y = 23.01 (1 -0.39 e -0.01389age ) 3 ) by Von Bertalanffy, y = 54.65 (1 -0.37 e -0.01577age ) and y = 58.91 (1 -0.36 e -0.01647age ) by Brody were fitted well to describe body weight, wither height and hip height of Kejobong goat, respectively.