Variants in GH, IGF1, and LEP genes associated with body traits in Santa Inês sheep

Growth hormone (GH), insulin-like growth factor 1 (IGF1), and leptin (LEP) can be candidate genes for association studies because they play vital roles in the metabolism process. Thus, this study aimed to identify variants in these genes associated with body traits in Santa Inês sheep. The following were recorded: body weight at 100 (BW100) and 240 days (BW240), average daily gain (ADG), withers (WH) and croup (CH) heights, body length (BL), thoracic (TG) and leg (LG) girths, thoracic (TW) and croup (CW) widths, body depth (BD), rib eye area (REA), fat thickness (FT), and carcass finishing score (CFS). Single-locus association analysis was performed with 11 variants in IGF1, 18 in LEP, and 16 in GH. Moreover, two haplotypes in IGF1 and one haplotype in LEP were evaluated in haplotype association analysis. The singlelocus analysis revealed 23 suggestive additive effects (p < 0.05), but no additive effect was found at the Bonferroni threshold. Haplotype association analysis revealed 19 additive effects, of which ten were at the Bonferroni threshold (p < 0.0074). In IGF1 gene, haplotype replacements were associated with ADG 20.51(7.37), CH 4.09(1.21), WH 3.52(1.20), BL 3.94(1.19), TG 3.88(1.30), TW 1.13(0.36), and LG 3.40(1.08); while in the LEP gene the haplotype replacement was associated with BW100 1.83(0.51), BD -2.51(0.56), and CFS -0.24(0.06). Therefore, there are haplotypes in IGF1 and LEP genes associated with body traits in Santa Inês sheep, which can be useful in marker-assisted selection.


Introduction
The Callipyge, Carwell, Double Muscling (Cockett et al., 2005), and Myostatin (Hickford et al., 2010) are major genes that affect growth and carcass traits in sheep. However, the variants in these genes do not explain total genetic variation related to either growth or carcass traits in sheep. Additionally, certain polymorphisms in these genes are typical of a number of specific sheep breeds (Cockett et al., 2005). Therefore, it is essential to study other candidate genes. In this context, the growth hormone (GH), insulin-like growth factor 1 (IGF1), and leptin (LEP) genes can be candidates for association studies with either growth or carcass traits in sheep. The transcripts of these genes play vital roles in the metabolism process (Tuersunjiang et al., 2016). Consequently, these genes can explain the genetic variance of several phenotypes in livestock.
The GH gene encodes the growth hormone which, through its specific receptor (Sahu et al., 2017), is the major regulator of IGF-I synthesis in the liver (Teran et al., 2016). IGF-1 is the main mediator of the effects of the growth hormone on muscle and bone tissues. Therefore, the GH/IGF is the major axis controlling the growth of animals (Yakar et al., 2018). A number of polymorphisms in GH were associated with growth (Abdelmoneim et al., 2017) and carcass traits (Gorlov et al., 2017) in sheep. In addition, associations between variants in the IGF1 gene with sheep growth traits were reported (Trukhachev et al., 2016).
Leptin has multiple functions, such as regulating appetite, body weight, maturation of the reproductive axis, neuroendocrine adaptations to fasting, and regulating glucose homeostasis (Caron et al., 2018). The vital role of leptin in controlling appetite can cause an impact on the variables directly associated with food intake, such as body weight and fat deposition. Thus, variants in the LEP gene were associated with growth (Hajihosseinlo et al., 2012), morphometric (Bakhtiar et al., 2017), and carcass traits (Barzehkar et al., 2009) in different sheep breeds.
Certain variants in the IGF1 gene were associated with internal carcass length, rib yield, and neck weight in Santa Inês sheep; while a variant in the GH gene was associated with weights and yields of primal cuts such as rib, loin, leg, and neck (Meira et al., 2019). Additive effects of polymorphisms in the LEP gene were also reported for weights and yields of these primal cuts, as well as weights and yields of the carcass (Meira et al., 2018). However, growth, morphometric, and in vivo ultrasound carcass traits were not evaluated in previous association studies on IGF1, GH, and LEP genes in Santa Inês sheep. Thus, this study aimed to identify polymorphisms in the GH, IGF, and LEP genes associated with body weight, morphometric traits, and carcass in vivo traits in Santa Inês sheep. Deoxyribonucleic Acid (DNA) was extracted from 192 Santa Inês lambs, all males, of which 74 were born in 2010, 15 in 2011, and 17 in 2012. The other 86 lambs were raised in 2014 in São Gonçalo dos Campos, BA, Brazil (12°25'58" S, 38°58'01" W, altitude of 233 m). All animals received diets formulated according to the National Research Council (NRC, 2007) to meet the nutritional requirements for lambs with estimated weight gains of 170 g d -1 .
The morphometric traits were measured in lambs at 240 days of age using a tape and a measuring stick. Withers (WH) (distance from the highest point of the thoracic vertebra to the ground) and the croup (CH) (distance from the coxal tuberosity to the ground) heights; body length (BL) (distance from supraglenoid tuberosity to sciatic tuberosity); thoracic (TW) (distance between the supraglenoid tuberosities) and croup (CW) (distance between the coxal tuberosities) widths; thoracic (TG) (contour of the thoracic cavity, adjacent to the shoulder blades) and leg (LG) (mid-thigh contour) girths; and the body depth (BD) (distance from the thoracic vertebrae to the sternum) were the morphometric traits evaluated.
The carcass finishing score (CFS) was accessed in vivo with values from 1 to 5, by a single recorder as follows: score 1 (very lean carcass), 2 (lean), 3 (medium), 4 (fat), and 5 (very fat). Additionally, ultrasonography was performed to obtain the rib eye area (REA) and fat thickness (FT) at 240 days of age. The REA and FT ultrasound images, between the 12 th and 13 th ribs on the left side of the lambs, were recorded. The length (A) and maximum depth (B) of the muscles were measured for estimating the REA using the equation A B 2 2 * * π . Furthermore, the average daily gain (g d -1 ) was calculated using the difference between the body weights obtained at 100 (weaning) and 240 days of age, divided by the number of days between these measurements. A descriptive statistical summary of the traits can be found in Table 1.

Genotyping
For the extraction of DNA a sample of 5 mL of blood was collected in vacutainer tubes containing Ethylenediaminetetraacetic acid (EDTA). The DNA extraction was performed using a salt precipitation method and proteinase K solutions following the Embrapa protocol (Oliveira et al., 2007). The primer design for amplification of the genes was determined by observing the available sequence in the National Center for Biotechnology Information (NCBI) database for LEP (GeneID:443534), IGF1 (GeneID:443318), and GH (GeneID:443329) genes of the sheep genome (version Oar_v4.0). The oligonucleotides were designed using the Primer 3 software package (http://bioinfo. ut.ee/primer3-0.4.0/). Net Primer was used to test the quality of the sequences. After selecting the forward and reverse primers, the Basic Local Alignment Search Tool (BLAST) was used for sequence alignment in NCBI (http://blast.ncbi.nlm. Nih.gov/Blast.cgi) and to confirm the similarity with the Ovis aries sequence. The forward (GTGCTGCTTTTGTGATTTCTTG) and reverse (GATAGAAGAGATGCGAGGAGGA) primers were used for the amplification of 4,550 bp of IGF1 between the positions 171037883 and 171113228. For the GH gene, the forward (GCTGCTGACACCTTCAAAGA) and reverse (TGACCCTCAGGTACGTCTCC) primers were used for the amplification of 1,194 bp between the positions 47485651 and 47487536. For the LEP gene, the forward (GGACCCCTGTACCGATTCCT) and reverse (CAAACTCAGGAGAGGGTGGA) primers were used to amplify 2,045 bp located between the positions 92501195 and 92503239.
Touchdown-PCR were performed on the three genes, using 15 µL of the reaction mixture, containing 0.3 mM of each of the primer, Taq polymerase, and 100 ng of the template DNA. The PCR conditions for LEP gene were as follows: initial denaturation of 98 °C for 5 min, followed by 20 cycles with denaturation at 98 °C for 10 s, annealing at 67 °C down to 57 °C, varying at -0.5 °C at each cycle for 30 s, and extension at 72 °C for 3 min. Immediately after the initial cycles, another 20 cycles were followed with denaturation at 98 °C for 10 s, followed by annealing at 57 °C for 30 s, extension at 72 °C for 3 min, and a final extension of 72 °C for 5 min. For the IGF1, an initial denaturation 98 °C/5 min, followed by ten denaturation cycles at 98 °C/10 s, annealing at 61 °C to 56 °C, reducing -0.5 °C at each cycle, for 30 s, and extension at 72 °C/4 min. Following another 30 cycles with denaturation temperature at 98 °C/10 s, annealing at 56 °C/30 s, and extension at 72 °C/4 min, ending with an extension at 72 °C/5 min. For GH amplification, the initial denaturation at 98 °C/5 min, followed by 20 denaturation cycles at 98 °C/10 s, annealing at 63 °C to 53 °C, reducing -0.5 °C at each cycle, for 30 s, and extension at 72 °C/1 min. A further 20 cycles with denaturation temperature at 98 °C/10 s were followed by annealing at 53 °C/30 s, and extension at 72 °C/1 min, ending with an extension at 72 °C/5 min. After PCR, the amplicons were purified with magnetic beads, and the recommended volume of AgencourtAMPure XP (Beckman Coulter, Brea, USA) was used to homogenize the beads to bind to the amplified products. Immediately after this step, the samples were purified with 70 % ethanol to remove the contaminants. The samples were quantified with Qubit® fluorometer (Life Technologies, Carlsbad, USA) and diluted to 0.2 ng µL -1 for library preparation. The Nextera® XT DNA sample preparation and the Nextera® XT index (Illumina, San Diego, USA) were used to prepare the library. All steps performed followed the Nextera XT protocol. Sequencing was performed on the MiSeq platform (Illumina, San Diego, USA) using the MiSeq Reagent Kit v2 (500 cycles).
The qualities of the reads were verified by the FastQC software program (https://dnacore.missouri. edu/PDF/FastQC_Manual.pdf). For the first data filtering, the SeqyClean software tool version 1.3.12 (Zhbannikov et al., 2013) was used, adopting a quality parameter of 24 (Phred score) for each base and a minimum length of 50 bp. Subsequently, the reads were aligned against the reference sheep genome deposited in the NCBI (version Oar_v4.0) by using the Bowtie2 program (Langmead and Salzberg, 2012). Finally, the functional annotation was performed using the variant effect predictor (VEP) for the online annotation of Ensembl in order to identify the locations of the mutations across different regions of the genome and the possible functional effects of the variants. For the nomenclature of the variants, we followed the recommendations of the Human Genome Variation Society (HGVs).

Haplotype and Hardy-Weinberg Equilibrium
The Hardy-Weinberg Equilibrium (HWE) was tested by comparing the predicted and observed heterozygosities. The predicted heterozygosity (PH) was obtained using the equation: PH = 2*(1 -MAF)* MAF, where MAF was the minor allelic frequency. The Haploview software (version 4.2) was used to test the HWE and to search for linkage disequilibrium (LD) blocks. The haplotype analysis revealed two LD blocks in IGF1 and one LD block in LEP, but no LD blocks were found in the GH gene.

Association analysis
Before the association analysis, all traits had been analyzed using the model: where y ijkl is the phenotypic value, µ the general average, F i the farm effect, Y j the year of birth effect, M k the month of birth effect, α ijkl (A) the covariate animal's age, and ε ij the error term. This analysis aims to identify possible record errors and test the assumptions of the analysis of variance (ANOVA). The MIXED procedure in SAS (Statistical Analysis System, version 9.3) was used in ANOVA. All the fixed effects tested were significant (p < 0.01). Therefore, it is important to include them in the association analysis models. Also, a principal component analysis to evaluate structuration was performed as reported by Price et al. (2006). The result (Figure 1) did not indicate structuration.
The single-locus association analysis was carried out using the software Qxpak 5 tool (Pérez-Enciso and Misztal, 2011), which performs a likelihood ratio test. The general model can be described as , where y is the vector containing the records of the traits, β the vector of solutions for the fixed effects, δ the vector of solutions for the genetic effects for any of the n QTL (quantitative trait loci) that affect the trait, and ε the vector of the residues. X and Z are the incidence matrices that correlate observations in y to the solutions in the vectors in β and δ k , respectively. The fixed effects included in the model were: (a) farm (2 levels), (b) year (4 levels), (c) month of birth (12 levels), and (d) the covariate age of the animal. The additive and dominance effects of the QTLs were also tested. The additive effect was calculated as the contrast between the genotypes (SS -RR), where the allelic variant (R) was found in the reference gene sequence, while (S) is the allelic variant found in Santa Inês sheep. The dominance effect was calculated as the contrast Only variants with MAF ≥ 2 % and HWE (p ≥ 0.05) were used in this analysis. For the single-locus analysis, forty-five variants (18 in LEP, 11 in IGF1, and 16 in GH) that showed HWE (p > 0.05) and MAF ≥ 1 % were employed (Table 2). Thus, the Bonferroni threshold for single-locus analysis was 0.0005.
The haplotype association analysis was performed as reported by Lake et al. (2003), and the haplo.glm subroutine of the haplo.stat package (version 1.7.7) was utilized in this analysis. Only haplotypes with a frequency greater than 4 % were tested. The haplotype analysis revealed two LD blocks in IGF1, one in LEP, and no LD blocks were found in the GH gene. In the haplotype association analysis, three LD blocks were used (Table 3). Therefore, the Bonferroni threshold was 0.0074.

Associations analysis -GH gene
Haplotype association analysis was not carried out for GH because no LD blocks were detected in this gene. Additionally, after Bonferroni correction, no association was found with the single-locus analysis approach in GH. However, the variant rs589527314, a substitution C/A found in the intron-2 of the GH gene, had a suggestive additive effect (p < 0.05) on BW100 (Table 4). The A allele was associated with a higher value of these traits, and the difference between the CC and AA genotypes was 3.2 kg. Associations analysis -IGF1 gene Single-locus analysis in IGF1 did not find an additive effect at the Bonferroni threshold. However, suggestive additive effects (p < 0.05) of the variant rs412470350 on WH, CH, TW, LG, and ADG were found (Table 4). This variant is a C >T substitution in the intron-1 of IGF1 gene. For all traits, the C allele was associated with a higher mean value. The differences  0.0088 LRT = likelihood ratio test; BW100 = Body weight at 100 days of age; WH and CH = withers and croup heigths; TW and CW = thoracic and croup widths; LG = leg girth; ADG = average daily gain; CFS = carcass finishing score; FT = fat thickness. Variants associated with body traits Sci. Agric. v.78, n.3, e20190216, 2021 between the CC and TT genotypes were: 1.88 cm (WH), 2.86 cm (CH), 1.15 cm (TW), 2.66 cm (LG), and 29.48 g d -1 (ADG). In addition, twelve haplotype replacements in LD block-1 of IGF1 were found (Table 5). The replacement of TCC by TCT was associated (p < 0.0075) with ADG, CH, WH, BL, TG, TW, and LG, where regression coefficients (β) and standard errors (SE) were 20.5079 ± 7.3741, 4.0859 ± 1.2121, 3.5227 ± 1.2002, 3.9393 ± 1.1920, 3.8814 ± 1.3036, 1.1302 ± 0.3577, and 3.3994 ± 1.0808, respectively. Moreover, another five suggestive (p < 0.05) haplotype replacements were found.

GH gene
The variant rs589527314 in intron-2 of the GH gene was associated with BW100 in Santa Inês sheep, which is supported by previous association studies. Variants in intron-2 were associated with body weight in Black Tibetan sheep (Han et al., 2016), and birth weight, body weight at 120 days, and ADG from birth to 120 days in Harri sheep (Abdelmoneim et al., 2017). Other regions of GH were also associated with body weight in sheep. The variant in intron-1 was associated with weaning weight and ADG in Nilagiri sheep (Cauveri et al., 2016). On exon 4, a PCR-SSCP associated with body weight in Makooei sheep was reported (Moradian et al., 2013), and one single nucleotide polymorphism was associated with body weight in Tibetan and Poll Dorset breeds (Jia et al., 2014). Moreover, a PCR-SSCP in exon 5 was associated with body weight in Baluchi sheep at six months of age (Tahmoorespur et al., 2011), while a PCR-RFLP was associated with body weight and ADG in Donggala and East Java breeds (Malewa et al., 2014) and Jambi Province sheep (Depison et al., 2017). These many associations between variants in GH and growth traits are supported by biological processes such as insulin secretion (Zhang et al., 2004), response to food (Berryman et al., 2004), and positive regulation of multicellular organism growth (Harvey, 2010), which depends on growth hormone activity.

IGF1 gene
The variant rs412470350 and haplotype replacements were associated with different body traits in Santa Inês sheep (Tables 4 and 5). The molecular function and biological process associated with the IGF1 gene supports these results. Yakar et al. (2002), demonstrated that circulating levels of IGF-1 directly regulate bone growth and density in mice, while Sun et al. (2014) found an association between IGF1 gene expression and body weight in Hu sheep. In addition, the IGF1gene has important molecular functions for animal growth, such as cell proliferation, maintenance of skeletal muscle satellite cells for regeneration, and an essential role in muscle hypertrophy (Philippou et al., 2007). Previous studies carried out with other sheep breeds also supported our results. A variant (rs430457475) in intron-1 was associated with WH and CH in Russian Merino breed (Trukhachev et al.,  LG = leg girth; BW100 and BW240 = Body weight at 100 and 240 days of age, respectively; REA = rib eye area; CFS = carcass finishing score; and BD = body depth. Variants associated with body traits Sci. Agric. v.78, n.3, e20190216, 2021 2016). Furthermore, association between a PCR-RFLP in the regulatory region of IGF1 with BL, TW, and body weight in Balami sheep, and WH in Yankasa sheep were also reported (Raji et al., 2017). A PCR-SSCP in exon-1 was associated with ADG and BL in Makui sheep , and ADG in both Baluchi (Tahmoorespur et al., 2009) and Makooei (Negahdary et al., 2013) sheep breeds.

LEP gene
Variants in the LEP gene of Santa Inês sheep were associated with several body traits in the current study. These results are supported by molecular functions such as hormonal activities (Wauters et al., 2000) that regulate vital metabolic processes; and the activity of growth factors that stimulate cells to grow or proliferate. Furthermore, the LEP gene is a component of several biological processes, especially those related to fat metabolism, such as the lipid metabolic process and beta-oxidation of fatty acids (Havel, 2004), which may also explain the associations observed here for FT and CFS. The leptin also plays a vital role in biological processes that are related to the negative regulation of appetite (Blundell et al., 2001), adult eating behavior and feeding behavior (Williamson et al., 2005), intestinal absorption (Yarandi et al., 2011), and bone growth (Upadhyay et al., 2015), which can explain the effects on growth and morphometric traits in Santa Inês sheep.
Previous studies also support our results with the LEP gene. In intron-2, associations are reported with REA in the Suffolk breed (Boucher et al., 2006) and with cold carcass weight, fat-tail, and total body fat in the Shal and Zelsheep breeds (Barzehkar et al., 2009). Additionally, the variants g.92501372G>A, rs398357543, rs421064645, and rs1135847360 in intron-2 of LEP gene were associated with several carcass traits recorded post-mortem in Santa Inês sheep, such as cold and hot carcass weights and yields, internal carcass length, leg and neck yields, neck weight, and CFS (Meira et al., 2018). Polymorphisms in other regions of the LEP gene, especially in exon 3, were also associated with economic traits. A PCR-SSCP in exon 3 of LEP gene was associated with breeding values for TG and rump length in Makooei sheep (Sadeghi et al., 2014), body weight and ADG in Makooei sheep (Hajihosseinlo et al., 2012), body weight at 90 days in Baluchi sheep (Tahmoorespur et al., 2010), and body weight in Kermani sheep with 3, 6, 9, and 12 months of age (Shojaei et al., 2010). Some variants in the exon 3 were associated with BL in Sanjabu sheep (Bakhtiar et al., 2017) feed conversion rate, residual feed intake and feed intake in crossbreed Awassi-Merino sheep (Jonas et al., 2016), and cold carcass weight in Santa Inês sheep (Quirino et al., 2016).
The present study found evidence of intronic variants in the IGF1, LEP, and GH genes associated with body traits in Santa Inês sheep. No previous studies identified small RNA transcripts encoded by the intronic regions associated with body traits in the current study. Moreover, the variants, here associated with body traits, are not in splice sites. Therefore, linkage disequilibrium between the polymorphism found here, and the causal variant is the main hypothesis that explains the additive effects found.