Genetic Polymorphisms of IGF1 and IGF1R Genes and Their Effects on Growth Traits in Hulun Buir Sheep

The identification of candidate genes and genetic variations associated with growth traits is important for sheep breeding. Insulin like growth factor 1 (IGF1) and insulin like growth factor 1 receptor (IGF1R) are well-accepted candidate genes that affect animal growth and development. The current study attempted to assess the association between IGF1 and IGF1R genetic polymorphisms and growth traits in Hulun Buir sheep. To achieve this goal, we first identified three and ten single nucleotide polymorphisms (SNPs) in exons of IGF1 and IGF1R in Hulun Buir sheep and then constructed six haplotypes of IGF1R based on linkage disequilibrium, respectively. Association studies were performed between SNPs and haplotypes of IGF1 and IGF1R with twelve growth traits in a population encompassing 229 Hulun Buir sheep using a general linear model. Our result indicated three SNPs in IGF1 were significantly associated with four growth traits (p < 0.05). In IGF1R, three SNPs and two haplotype blocks were significantly associated with twelve growth traits (p < 0.05). The combined haplotype H5H5 and H5H6 in IGF1R showed the strong association with 12 superior growth traits in Hulun Buir sheep (p < 0.05). In conclusion, we identified SNPs and haplotype combinations associated with the growth traits, which provided genetic resources for marker-assisted selection (MAS) in Hulun Buir sheep breeding.


Introduction
Growth traits are among the most important economic attributes in sheep breeding and are of great concern to breeding experts. Growth traits, including body weight, average daily gain and body size greatly influence meat productivity, which influences production and profitability in the mutton sheep industry [1]. Studies have revealed that many candidate genes are related to growth traits, among which IGF1 and IGF1R genes are wellaccepted candidate genes that affect growth and production performance in livestock [2,3]. Insulin-like growth factor 1 (IGF1) is an endocrine growth factor involved in normal growth and development [4][5][6], fetal development and metabolism [7,8]. Insulin-like growth factor 1 receptor (IGF1R) is encoded by the IGF1R gene and is a receptor tyrosine kinase that mediates the actions of IGF1 [9,10].
Significant associations were identified between single nucleotide polymorphisms (SNPs) of the two genes and growth performance in diverse farm animals, including cattle [11][12][13], buffaloes [14], pigs [15,16] and goats [17][18][19][20]. In sheep, it has been reported that SNPs of IGF1 and IGF1R are related to meat production and growth [21][22][23]. Using the PCR restriction fragment length polymorphism (PCR-RFLP) method, Grochowska et al. found a highly significant effect of SNPs in the 5' untranslated (5' UTR) region of IGF1 on carcass traits and meat compositions in local sheep breeds in Poland Merino sheep [24]. Negahdary et al. found a significant effect of the 5' UTR region of the IGFI gene on birth weight, weaning weight, 6-month weight, average daily gain from birth to weaning and average daily gain from 6 to 9 months in Makooei sheep [25]. A mutation in intron 12 of the IGF1R gene was significantly associated with body weight and growth rate in Pomeranian Coarse wool ewes [26]. Later, associations were found between SNP in exon 3 of IGF1R and daily gain in the early developmental stage of Colored Polish Merino sheep [23]. The discovery of associations between genetic polymorphisms and growth traits provides useful information for the genetic improvement in sheep breeding. SNPs in the exon of genes are important because they may cause potentially functional variations, which lead to phenotypic changes in livestock. Since most identified SNPs in IGF1 and IGF1R were located in the 5' flanking regions, we paid particular attention to genetic variations in the exons of the two genes.
Hulun Buir sheep are one of the representative indigenous sheep breeds in northern China, characterized by their high-grade meat quality and outstanding resistance to stress, such as cold and roughage. A lack of advanced breeding methods leads to poor growth performance compared to commercial breeds. To conduct the genetic improvement and breeding in Hulun Buir sheep, several works have been performed to identify genetic variations that were associated with economic traits in Hulun Buir sheep. It has been shown that the somatostatin receptor 1 (SSTR1) gene harbors two SNPs that were remarkably associated with growth traits of Hulun Buir sheep [27]. Based on Genome-wide association studies (GWAS), six SNP loci from 526,225 autosomal markers were greatly associated with carcass traits and chest girth [28]. Candidate genes and SNPs have been reported to be associated with fat deposition and fat metabolism [29,30]. However, no systematic investigations have been reported on the association between genetic polymorphism and the early growth traits in Hulun Buir sheep.
To improve the growth performance of Hulun Buir sheep, we investigated the genetic polymorphisms of IGF1 and IGF1R and their associations with twelve growth traits. By scanning exons of IGF1 and IGF1R, we identified thirteen SNPs in the IGF1 and IGF1R genes and two haplotype blocks involving six haplotypes in 229 Hulun Buir sheep. Among these SNPs and haplotype blocks, six SNPs and two haplotype blocks were remarkably associated with growth traits in Hulun Buir sheep. Notably, we identified two combined haplotypes that demonstrated a strong association with twelve greater phenotypic traits. Conclusively, our study provided useful information and laid the foundation for the genetic breeding of Hulun Buir sheep.

Animals and Data Collection
In total, 229 Hulun Buir lambs (male = 106, female = 123), which were born in March 2019 on Hulun Buir sheep farms (Hulun Buir city, Inner Mongolia, China), were investigated. The animals were grazed in identical conditions. The birth weight (BW), weaning weight adjusted at 4-month-old (WW) and body weight at 9-month-old (NBW) were recorded. Meanwhile, average daily gains (ADG) during birth to weaning, weaning to 9-month-old and birth to 9-month-old periods were calculated. Body height (BH), body length (BL) and chest girth (CG) were measured at weaning and at 9 months of age, respectively. Approximately 1 cm 3 marginal ear tissues were collected and preserved in 95% ethanol. Genomic DNA of Hulun Buir sheep was extracted using a Tiangen DNA extraction kit (Tiangen Biotech Co., Ltd, Beijing, China) and stored at −20 • C for PCR amplification.

SNP Identification and Genotyping
Four and twenty-one pairs of primers were designed for all exons of IGF1 and IGF1R genes based on the published mRNA sequences (Gene ID: 443318, GenBank No. NM_001009774 (IGF1), Gene ID: 443515, GenBank No. XM_027957015 (IGF1R)), using Primer Premier 5.0 (Premier Biosoft, Palo Alto, Santa Clara, CA, USA), respectively. The primer information is listed in Supplementary materials Table S1. The PCR contained 100 ng template DNA, 10 pM of each primer, 3.5 µL 10 × PCR buffer, 2.5 mM dNTP, 1 U of Taq DNA polymerase (Takara Biotechnology Co., Ltd., Beijing, China) and double-distilled water (ddH 2 O), to make up a volume of 35 µL. PCR was performed in Thermocycler System (ABI 9700, Applied Biosystems, Waltham, MA, USA) with the following reaction procedure: predenaturation at 94 • C for 5 min, followed by 35 cycles at 94 • C for 30 s, 53-60 • C for 30 s and 72 • C for 40 s, with a final extension at 72 • C for 10 min. PCR products were separated by gel electrophoresis (1.5% agarose), purified using magnetic beads (Agencourt AMPure XP, Beckman Coulter, Krefeld, Germany) and sequenced in an Agilent 3730 sequencer (Agilent Technologies, Santa Clara, CA, USA). The sequencing results were aligned to published sheep IGF1 and IGF1R genes using Chromas 2.0 and SeqMan (DNASTAR software, version 7.1) to identify potential SNPs.

Population Genetics of IGF1 and IGF1R Genes
Genotypic and allelic frequencies were estimated with the direct counting method. Hardy-Weinberg equilibrium (HWE), observed heterozygosity (Ho), expected heterozygosity (He) and effective allele numbers (Ne) were analyzed according to the genotype frequencies of SNPs [31]. Cervus (version 3.0) was used to calculate the polymorphic information content (PIC) of each mutation site [32].

Linkage Disequilibrium Analysis and Haplotype Construction
The extent of linkage disequilibrium (LD) between each pair of SNPs in IGF1 and IGF1R was analyzed according to the value of r 2 using Haploview software (version 4.2) [33]. Haplotype blocks with strong LD of SNPs (r 2 > 0.33) were defined based on the confidence intervals methods [34].

Statistical Analyses
SAS software (version 13.0, SAS Institute) was applied for statistical analyses, and the results were expressed as the mean ± SE (standard error). The associations were carried out between the genotypes and individual growth traits using general linear model (GLM): where Y ij is a growth trait measured on an individual animal (BW, WW, NBW, ADG, BH, BL and CG); µ is the mean value of overall; G i is the fixed effect of genotypes of the population (i = 3 levels, except rs600896367 of IGF1 gene and c.244C>T, rs162159917, rs601806812 and rs193644211 of IGF1R gene where i = 2 levels); S j is the fixed effect of sex (j = 2 levels); G i × S j is the interaction effect between sex and genotypes; if the difference of interaction effect between the sex and genotypes is not significant, the general linear model should be reduced as: The association analysis between haplotype combinations and individual growth traits was analyzed by the following GLM: where Y ij is a growth trait measured on an individual animal (BW, WW, NBW, ADG, BH, BL and CG); µ is the mean value of overall; H i is the fixed effect of haplotype combinations of the population (i = 5 levels); S j is the fixed effect of sex (j = 2 levels); H i × S j is the interaction effect between sex and haplotype combinations; if the difference of interaction effects between the sex and haplotype combinations is not significant, the general linear model should be reduced as: ε is the random error in the above models. Tukey's test and Bonferroni corrections were performed for multiple pairwise comparisons between genotypes or haplotype combinations based on SNPs. The p value of 0.05 was defined as statistical significance.

SNP Detection of IGF1 and IGF1R Genes in Hulun Buir Sheep
We detected three SNPs in the IGF1 gene and ten SNPs in the IGF1R gene in 229 Hulun Buir sheep ( Figure 1, Figure 2, Table 1). All of the detected SNPs were transition mutations except for SNP13 (transversion mutation) in exon 19 of IGF1R. SNP4 in IGF1R was a nonsynonymous mutation resulting in a substitution of Cys for Arg in amino acid sequence, and the rest were synonymous mutations. By searching in the dbSNP database of NCBI, we found that SNP4 and SNP8 in IGF1R were two novel single-nucleotide mutations in sheep and will be uploaded to the SNP data bank (Table 1).
Hulun Buir sheep ( Figure 1, Figure 2, Table 1). All of the detected SNPs were transition mutations except for SNP13 (transversion mutation) in exon 19 of IGF1R. SNP4 in IGF1R was a nonsynonymous mutation resulting in a substitution of Cys for Arg in amino acid sequence, and the rest were synonymous mutations. By searching in the dbSNP database of NCBI, we found that SNP4 and SNP8 in IGF1R were two novel single-nucleotide mutations in sheep and will be uploaded to the SNP data bank (Table 1).
SNP13 rs161167008 exon19 C G Gly synonymous  the arrow indicates the T-C mutation site. SNP3: c.495G>A (rs400398060); the arrow denotes to the G-A mutation site.

Genotyping, Genotypic and Allelic Frequencies
Among all the SNPs, the wild types were dominant alleles compared with the mutants ( Table 2). Genotyping results showed that SNP1 in IGF1 as well as SNP4, SNP5, SNP9 and SNP12 in IGF1R displayed two genotypes: wild-type homozygotes and mutant heterozygotes, and the remaining eight SNPs showed three different genotypes: wild-type homozygotes, mutant heterozygotes and mutant homozygotes ( Table 2). In SNP6-8, heterozygotes showed the highest genotype frequencies compared with wild-type and mutant homozygous. In the remaining 10 SNPs, the wild-type homozygotes had the highest genotype frequencies compared with mutant heterozygotes and homozygotes ( Table 2).

Genetic Diversity and Hardy-Weinberg Equilibrium
The allelic frequencies of all 13 SNPs obey the HWE law (p > 0.05). The Ne values of SNP2 in IGF1 and SNP6-SNP8 in IGF1R were close to 2. The PIC value showed that the five SNP loci (SNP2, SNP3, SNP6-SNP8) exhibited low polymorphism (PIC < 0.25), while the remaining eight SNPs showed moderate polymorphism in the Hulun Buir sheep population (0.25 < PIC < 0.5) ( Table 2).

Effects of Genotypes on Growth Traits
Association analysis was performed between genotypes of the SNPs and growth traits on 229 Hulun Buir sheep. The statistical results were listed in Supplementary Tables S2-S5.

Effects of SNP Genotypes in IGF1 on Growth Traits
The GA genotype of SNP1 had significantly greater WCG and NBL than the GG genotype (p < 0.05). At the SNP2 locus, the higher NCG was observed in TC genotype than that in the CC genotype but not in the TT genotype (p < 0.05). The GG and GA genotypes of SNP3 were significantly associated with greater 4-9 ADG than the AA genotype (p < 0.05, Figure 3).  18   The GA genotype of SNP1 had significantly greater WCG and NBL than the GG 24 genotype (p < 0.05). At the SNP2 locus, the higher NCG was observed in TC genotype than 25 that in the CC genotype but not in the TT genotype (p < 0.05). The GG and GA genotypes 26 of SNP3 were significantly associated with greater 4-9 ADG than the AA genotype (p < 27 0.05, Figure 3).  33 SNP3 genotypes of IGF1 gene; 4-9 ADG = average daily gain from 4 to 9-months-old. Different 34 letters (small letters: p < 0.05) above the column indicate significant differences among the different 35 genotypes. 36 3.3.2. Effects of SNP Genotypes in IGF1R on Growth Traits 37 The mutant homozygotes (CC) of SNP6 had significantly longer NBL than those 38 individuals with the TC genotype (p < 0.05, Figure 4A). Significant differences (p < 0.05) 39

Effects of SNP Genotypes in IGF1R on Growth Traits
The mutant homozygotes (CC) of SNP6 had significantly longer NBL than those individuals with the TC genotype (p < 0.05, Figure 4A). Significant differences (p < 0.05) and extremely significant differences (p < 0.01) were found between genotypes of the SNP8 locus with the 11 growth traits out of 4-9 ADG (Figure 4B,C). The genotypes containing the wild-type allele had better phenotypic values than mutant homozygotes. At the SNP13 locus, the individuals with the CC genotype had greater NBW, 0-9 ADG, WCG, NBH and NCG than those with the GG genotype (p < 0.05); the CC and CG genotypes were associated with significantly longer NBL than the GG genotype (p < 0.05, Figure 4D,E). No significant effects were detected among the remaining seven SNP loci and early growth traits of Hulun Buir sheep (p > 0.05).

Linkage Disequilibrium and Haplotype Analysis
A strong linkage disequilibrium (r 2 > 0.33) was observed among SNP5, SNP9 and SNP11, and between SNP6 and SNP7, as well as SNP8 and SNP9 loci in the IGF1R gene ( Figure 5). In particular, SNP6 to SNP9 loci formed two haplotype blocks. The first haplotype block was composed of SNPs 6 and 7, including three common haplotypes. The haplotypes H1 (TC), H2 (CT) and H3 (CC) occurred at frequencies of 0.537, 0.389 and 0.074, respectively, and five haplotype combinations were generated ( Table 3). The second haplotype block was composed of SNP8 and SNP9, including three common haplotypes. The haplotypes H4 (CG), H5 (TG) and H6 (CA) occurred at frequencies of 0.321, 0.581 and Genes 2022, 13, 666 8 of 14 0.098, respectively, and generated five haplotype combinations (Table 4). We did not detect the linkage disequilibrium among three SNP loci (r 2 < 0.33) in the IGF1 gene ( Figure 6).   Figure 5. Linkage disequilibrium plot (r 2 ) and haplotype blocks for SNPs of the IGF1R gene in Hulun Buir sheep. The values in boxes are pairwise SNP correlations (r 2 ); dark red boxes indicate strong LD (r 2 > 0.33) and light red boxes without numbers represent very weak LD (r 2 < 0.001).

Discussion
The growth of the animal was subject to growth hormone (GH)-IGF1 somatrotropic axis, in which GH acts as a major regulator for development, growth and anabolic processes. IGF1 modulates the biological actions of GH by binding to its receptor (IGF1R) [35]. IGF system includes IGF ligands and their receptors, which influences glycogenesis, glucogenesis and protein synthesis through the regulation of downstream gene expression and signaling pathways [36]. Among IGF ligands and receptors, IGF1 and IGF1R proteins are crucial regulators of cell growth and metabolism [37,38]. Genetic variation may have an impact on the phenotypic characteristics of animals by influencing The comparison of body weight traits for the haplotype combinations (block 2) of IGF1R gene in Hulun Buir sheep; BW = birth weight; WW = weaning weight (4-month-old); NBW = body weight at 9-months-old; 0-4 ADG = average daily gain from birth to 4-months-old; 4-9 ADG = average daily gain from 4 to 9-months-old; 0-9 ADG = average daily gain from birth to 9-months-old. (C) Association analyses for the haplotype combinations (block 2) of IGF1R gene with body size traits in Hulun Buir sheep; WBH = body height at 4 months of age; WBL = body length at 4-months-old; WCG = chest girth at weaning (4-month-old); NBH = body height at 9-months-old; NBL = body length at 9-months-old; NCG = chest girth at 9-months-old. Different letters (small letters: p < 0.05; capital letters: p < 0.01) above the column indicate significant differences among the different haplotype combinations.

Discussion
The growth of the animal was subject to growth hormone (GH)-IGF1 somatrotropic axis, in which GH acts as a major regulator for development, growth and anabolic processes. IGF1 modulates the biological actions of GH by binding to its receptor (IGF1R) [35]. IGF system includes IGF ligands and their receptors, which influences glycogenesis, glucogenesis and protein synthesis through the regulation of downstream gene expression and signaling pathways [36]. Among IGF ligands and receptors, IGF1 and IGF1R proteins are crucial regulators of cell growth and metabolism [37,38]. Genetic variation may have an impact on the phenotypic characteristics of animals by influencing the expression and function of the genes [39,40]. Therefore, we inferred that the genetic variation in IGF1 and IGF1R may also influence the growth traits of sheep.
In the present study, we discovered genetic polymorphisms of the IGF1 and IGF1R genes and evaluated their effects on growth traits in Hulun Buir sheep. Our results indicated that IGF1 and IGF1R exhibited low to medium genetic diversity, and some of the genetic variations exhibited a significant association with the growth performance in Hulun Buir sheep. This observation provided SNP marker information, which has potential feasibility for MAS in Hulun Buir sheep breeding schemes.
The Hardy-Weinberg equilibrium of all 13 SNPs indicated the absence of artificial selection of Hulun Buir sheep [41]. In the current study, two novel single-nucleotide polymorphisms were identified, including a nonsynonymous mutation of SNP4. A growing body of evidence has shown that the synonymous mutations could influence phenotypic performance by influencing gene expression through the regulation of mRNA stability and protein expression [42][43][44][45]. Maria et al. reported that the synonymous mutation rs159876393 SNP1 of IGF1 was associated with milk protein and casein contents in Sarda sheep [46]. A synonymous mutation SNP2 (rs159876393) in exon 2 of IGF1 was associated with variations in carcass traits of New Zealand Romney Sheep, including carcass weight, backfat thickness and the lean meat percentage [47]. Consistent with previous reports on other sheep breeds, we also identified a strong association of SNP1 and SNP2 with growth traits in Hulun Buir sheep, which indicated that SNP1 and SNP2 of the IGF1 gene might be related to multiple traits in sheep. A remarkable association was found between SNP3 (rs400398060) of the IGF1 gene and average daily gain from 4-9 months of age (4-9 ADG) in the present study. This mutant site was also detected in Egyptian Barki sheep and was not correlated with growth traits, indicating that its association might be dependent on the genetic backgrounds of sheep breeds [48]. Few studies reported the association between genetic polymorphisms of the IGF1R gene and growth traits in sheep. A significant correlation was detected between average daily gain and an SNP of the IGF1R gene in local sheep breed in Poland Merino sheep [23]. The present study reported 10 SNPs in the IGF1R gene, and SNP6, SNP8 and SNP13 were significantly associated with growth traits in Hulun Buir sheep. In addition, the sheep with homozygous wild genotype TT of SNP8 and CC of SNP13 had superior growth traits than those with homozygous mutant genotypes CC and GG, suggesting that they could serve as the predominant genotypes.
Generally, linked SNP loci are of much concern because of the existence of substantial LD between causal SNPs [49]. Haplotype combinations involving multiple linked SNP loci may provide more precise information than single SNP markers for association analysis [50][51][52]. In this study, the strong LD suggested that these alleles were tightly linked; thus, we carried out an association analysis between the haplotypes and growth traits. The association and multiple comparison analyses demonstrated that the H5H5 (TGTG) haplotype combination with wild-type alleles was the dominant haplotype. This was consistent with the result that the wild-type allele T of SNP8 was related to better growth traits. Additionally, SNP6-SNP9 formed two haplotype blocks, which displayed a remarkably significant effect on growth traits. Based on the results above, we inferred that the four SNPs did not act independently [53], and SNP6 and SNP8 of IGF1R may be causal mutations that affect phenotypic traits [54].

Conclusions
Conclusively, our analysis showed that SNP1, SNP2 and SNP3 of the IGF1 gene, SNP6, SNP8 and SNP13, as well as haplotype block 1 and haplotype block 2 of the IGF1R gene can be used as candidate markers for early growth traits in MAS of Hulun Buir sheep. Further studies will be conducted to investigate the effects of these SNPs on other economic traits in Hulun Buir sheep. The wild-type alleles of SNP8, haplotype combinations H5H5 (TGTG) and H5H6 (TGCA) in the IGF1R gene showed superior growth traits during the early stage.
Overall, our study provided important genetic variations, which could serve as potential markers for growth trait selection in Hulun Buir sheep.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13040666/s1, Table S1: Primer information of IGF1 and IGF1R in Hulun Buir sheep; Table S2: Associations for the SNPs of IGF1 gene with body weight traits and ADG traits in Hulun Buir sheep; Table S3: Associations for the SNPs of IGF1 gene with body size traits in Hulun Buir sheep; Table S4: Associations for the SNPs of IGF1R gene with body weight traits and ADG traits in Hulun Buir sheep; Table S5: Associations for the SNPs of IGF1R gene with body size traits in Hulun Buir sheep; Table S6: Associations for the haplotype combinations (block 1) of IGF1R gene with body weight traits and ADG traits in Hulun Buir sheep; Table S7: Associations for the haplotype combinations (block 1) of IGF1R gene with body size traits in Hulun Buir sheep; Table  S8: Associations for the haplotype combinations (block 2) of IGF1R gene with body weight traits and ADG traits in Hulun Buir sheep; Table S9: Associations for the haplotype combinations (block 2) of IGF1R gene with body size traits in Hulun Buir sheep.