Identification of Genes Targeted in South African Merino and Afrino Sheep Populations Under Long-Term Selection for Reproduction and Body Weight

Abstract

Background: Reproductive performance and body weight are of the utmost economic importance in determining the efficiency of sheep production. Simultaneous selection for increased reproductive performance and early growth traits is a common strategy in many flocks, but ambiguous results regarding the relationship between reproduction and body weight have been reported. The objective of this study was to perform a genome-wide association study (GWAS) in two South African Merino flocks and an Afrino sheep flock that were selected for both reproduction and body weight over decades. The GWAS aimed to identify SNPs associated with genes affecting the traits number of lambs born (NLB), number of lambs weaned (NLW), total weight of lamb weaned (TWW) and body weight (BW) and thus to ascertain which genes were targeted through directional selection.

Results: In the GWAS, 16 SNP markers associated with reproductive traits were identified among the three populations, while 15 SNPs were associated with body weight. These SNPs were linked respectively to 26 and 21 documented genes in the sheep genome. Most of these genes were previously associated in literature with reproduction related, as well as with growth related traits in various farm animal species. This study, supported by results from previous studies performed on sheep and cattle, identified the following genes that warrant further investigation as to their functions and processes relating to growth and reproduction in sheep: MAP7D1, TRAPPC3, THRAP3, TRMP8, SPP2, HDAC9, ZFHX3, SIX6, C14orf39, TAF4B, TRSP1 EYA2, RBMS3,  STL38L, BSPH1, LIG1, CABP5 and ELSPBP1.

Conclusions: Long-term selection in the flocks for both body weight and reproductive traits, and especially on the composite trait TWW, have favoured genes with pleiotropic effects influencing both groups of traits. SNPs associated with these pleiotropic genes were detected in the association analyses for the various traits.

Background

In most sheep breeding enterprises, reproductive performance and body weight are two main groups of traits that determine the efficiency of production. Apart from its direct impact (number of marketable lambs), a high reproductive rate also contributes to higher selection intensity. Furthermore, body weight is an important trait considered during selection of both replacement ewes and rams in many wool and mutton sheep production systems. Generally, the aim is to increase reproductive performance and early body weight and growth, while maintaining mature weight to limit maintenance requirements of the ewe flock under extensive conditions [1] (Herselman & Olivier, 2010). Consequently, simultaneous selection for both reproduction and body weight traits have to be performed. Although the genetic relationship between reproduction and body weight in most sheep flocks is positive [2–4] (Snyman et al., 1998; Safari et al., 2005; Olivier, 2014), negative relationships have also been reported [3] (Safari et al., 2005). Within flocks, variation among animals regarding the relationship between reproduction and body weight also exists. Various combinations of positive and negative breeding values for body weight and reproduction are present among animals in a specific flock, which may present difficulties in selecting breeding sires and dams.

One way of incorporating early growth and reproduction in a single selection criterion is through total weight of lamb weaned per ewe (TWW). TWW per ewe lifetime includes ewe fertility, litter size, lamb survival rate and direct growth performance of lambs until weaning and can be used as a biological index to measure reproduction potential in ewes [5, 6] (Zishiri et al., 2013; Matebesi-Ranthimo et al., 2017). TWW was used as reproductive selection criterion, together with body weight, in several governmental experimental sheep flocks in South Africa over many years. The data collected on these flocks will provide insight into the genomic consequences of simultaneous long-term selection for body weight and TWW as reproductive criterion.

Genome-wide association study (GWAS) is an important tool for the identification of candidate genes and molecular variants associated with economically important traits in farm animals. The identification of genes associated with reproduction and body weight in flocks subjected to simultaneous selection for these traits for many years, could provide information to elucidate the underlying mechanisms influencing these traits and the relationship between these traits. During the past decade, the number of GWAS in sheep has increased and various candidate genes have been identified for various body weight, growth and reproductive traits [7, 8] (Xu & Li, 2017; Gebreselassie et al., 2020). The identified genes varied between studies and breeds, and were identified across the entire genome.

The objective of this study was to perform a GWAS in two South African Merino flocks and an Afrino sheep flock that were selected for both reproduction and body weight for many years. The aim of the GWAS was to identify significant SNPs associated with genes affecting reproduction and body weight and thus to determine which genes were targeted when selection was based on body weight and the composite trait total weight of lamb weaned.

Results

Breeding value statistics

The average, minimum and maximum estimated breeding values (EBV) for each trait used in the analyses for the three populations are summarised in Table 1. There was a wide range in EBV among the animals within each flock in all the traits included in the GWAS.

Table 1 The average, minimum and maximum estimated breeding values (EBV) for each trait per population

Genetic relatedness within and between populations

The relatedness within and between the three populations were investigated and are illustrated in Figure 1. As the three populations clustered according to geographical region and showed genetic differentiation, each was included as a separate population in the GWAS analyses. The relatedness between the two Merino populations could be explained by the use of certain rams as sires in both populations.

Genome-wide association study

Manhattan plots of the four traits for the different populations are illustrated in Figures 2 to 5. Several suggestive SNPs were identified in all the populations for all the traits. These, as well as the SNP effects are summarised in Table 2. For number of lambs born, number of lambs weaned, total weight of lamb weaned and body weight, 10, 16, 9 and 14 SNPs were respectively identified across the three populations.

Table 2 Suggestive SNPs associated with reproduction and body weight in the three populations

The Q-Q plots for the various analyses are also illustrated in Figures 2 to 5. These plots for the Afrino body weight (Figure 2) and Cradock Merino NLB (Figure 3) indicate the largest deviations from the distribution under the null hypothesis, suggesting the strongest associations for these traits. Moderate deviations were observed for NLW in the Grootfontein and Cradock Merino populations (Figure 4) and for TWW in the Grootfontein Merino population (Figure 5).

The genes associated with the suggestive SNPs in the three populations are summarised in Table 3. Genes within 50 kb upstream or downstream of the SNP were included for the purpose of this study. Some of the suggestive SNPs were associated with more than one gene. MAP7D1, TRAPPC3 and THRAP3 are all located within 50 kb up- or downstream from SNP OAR1_10554666.1, while SPP2 and TRPM8 were both associated with SNP s20120.1 on OAR1 in the Cradock Merino flock. In the Afrino flock, RPB1 and RPB2 were located close to SNP OAR1_267463962.1. Two SNPs on OAR14 were also associated with more than one gene. ZNRF1, LDHD and ZFP1 were associated with SNP OAR14_2287469.1 and BSPH1, LIG1, CABP5 and ELSPBP1 with SNP OAR14_56900862.1.

Nine genes were associated with reproduction in the Cradock population. Two of these genes (GRIK3 and HDAC9) were associated with all three reproductive traits, while MAP7D1, TRAPPC3 and THRAP3 were associated with both NLB and NLW. Eight and nine genes were respectively associated with the reproductive traits in the Grootfontein and Afrino populations. The SIX6 and C14orf39 genes in the Afrino population were associated with all three reproductive traits, and C1orf68 was associated with both NLW and TWW. This illustrates the inter-relationships among these reproductive traits on a genomic basis. Two, six and 12 genes were associated with BW in the Cradock, Grootfontein and Afrino populations respectively.

Table 3 Genes located within 50 kb of the SNPs associated with reproduction and body weight in the three populations

Gene ontology (GO) categories identified for genes associated with body weight and reproduction in the Afrino, Grootfontein and Cradock populations are presented in Supplementary Tables 1 to 3 respectively. The gene groups identified in the various populations were classified into several main biological processes and molecular functions.

In Table 4 (Afrino population), the genes ENSOARG00000006463, C1orf68, C14orf39, MRPL2, CUL7 and ELSPBP1 were not included in any category, while the following 15 genes were categorised into the three GO categories: LIG1, SIX6, TAF4B, ZFP1, BSPH1, IL1RN, CABP5, KLC4, LDHD, PKP2, PTK7, RBP1, RBP2, STK38L, ZNRF1. Main biological processes identified were “Metabolic processes” and “Developmental processes”, while “Binding” was the only molecular function identified. Genes in this population that were categorised under biological processes were mostly associated with BW and TWW. All traits were associated with genes in the molecular function and cellular components categories.

Thirteen of the 14 genes identified in the Grootfontein population (EYA2, RBMS2, PSAT1, TRPS1, ZFHX3, ZNF704, PIK3C2A, XIRP2, CAP1, CEP78, ENKUR, ITFG1, THNSL1) were categorised into the three GO categories (Table 5). Only the PPT1 gene was not categorised. Main biological processes identified in the Grootfontein gene group were “Biosynthetic processes”, “Transcription”, “Metabolic processes” and “Gene expression”. All genes were categorised under “Binding” as molecular function. Similar to the Afrino population, genes categorised under biological processes were mostly associated with BW and TWW. The same applies for the “Binding” molecular function, with the addition of NLW. All traits were associated with genes under the cellular components category.

As summarised in Table 6, in the Cradock population, the following six of the 11 identified genes were categorised into the three GO categories, namely GRIK3, SLC4A10, THRAP3, TRPM8, HDAC9, ZNF496. The genes MAP7D1, TRAPPC3, SPP2, IER5L and ENSOARG00000005343 were not included in any category. Two biological processes, namely “Homeostasis” and “Transcription” and two molecular functions, “Transporter activity” and “Channel activity” were identified. All three reproductive traits were associated with genes categorised in all three GO categories, while BW was associated with genes in the “Transcription” biological process and cellular components category.

Table 4 Gene ontology categories for genes associated with body weight and reproduction in the Afrino population

^a BW = body weight; NLB = number of lambs born; NLW = number of lambs weaned; TWW = total weight of lamb weaned

Table 5 Gene ontology categories for genes associated with body weight and reproduction in the Grootfontein Merino population

^a BW = body weight; NLB = number of lambs born; NLW = number of lambs weaned; TWW = total weight of lamb weaned

Table 6 Gene ontology categories for genes associated with body weight and reproduction in the Cradock Merino population

^a BW = body weight; NLB = number of lambs born; NLW = number of lambs weaned; TWW = total weight of lamb weaned

Discussion

Extensive reviews on identified candidate genes associated with different growth and reproduction traits in various sheep breeds published by [7] Xu & Li (2017) and [8] Gebreselassie et al. (2020) confirmed that these genes are distributed throughout the sheep genome. There was limited agreement between the previously reported studies in either the identified genes or the broader genomic regions where the genes were located. This could be ascribed to various factors, including breed differences, long-term selection practices followed in the breeds or flocks, sample sizes, models of analyses applied, levels of significance for identification of significant and suggestive SNPs and distance of significant or suggestive SNPs from the associated gene. For example, [9] Al-mamun et al. (2015) identified 39 SNPs associated with body weight in Australian Merino sheep, with a region on OAR6 containing 13 significant SNPs. Several SNPs related to growth and carcass traits in Scottish Blackface lambs were also reported on OAR6 ([10] Matika et al., 2016). However, no significant or suggestive SNPs were found for body weight or growth traits on OAR6 by [11] Zhang et al. (2013) or [12] Almasi et al. (2020), nor in the populations in the current study.

Many of the genes associated with reproduction in the current study were previously associated in literature with reproductive related traits, as well as with growth related traits in different farm animal species. Such genes associated with reproductive traits in the current study are MAP7D1, TRAPPC3, THRAP3, SPP2 and HDAC9 in the Cradock Merino population, ZFHX3 and PIK3C2A in the Grootfontein Merino population and SIX6 in the Afrino population.

Previous literature linked many of the genes associated with body weight in the current study to reproduction traits in farm animal species. Some of these genes were associated with both growth and reproductive related traits in literature. In the case of the Afrino population, all the genes associated with body weight in the current study (except for LDHD), were previously related to reproductive traits (RBP1, RBP2, PKP2, STK38L, IL1RN, ZNRF1, ZFP1, BSPH1, LIG1, CABP5, ELSPBP1). The genes TRPS1, ENKUR, EYA2 and RBMS3 in the Grootfontein Merino population and ZNF496 in the Cradock Merino population were also previously associated with reproductive traits.

Long-term selection in the flocks for both body weight and reproductive traits, especially the composite trait TWW, could have favoured genes with pleiotropic effects influencing both these traits, with the result that SNPs associated with such pleiotropic genes would be detected in the GWAS for the different traits. Furthermore, genes associated with BW and TWW in the Afrino and Grootfontein Merino populations featured in the same GO biological processes identified for these two populations. This confirms that the same underlying physiological processes are involved in these traits.

Previously published literature linked genes associated with either reproduction or body weight in this study to various reproductive processes from oocytogenesis through conception, implantation and pregnancy to milk production. Genes differentially expressed in or associated with bovine follicles or oocytes were MAP7D1 and ZNF496 [13] (Dickinson, 2016), TRAPPC3 [14] (Donnison & Pfeffer, 2004), RBP1 [15] (Mamo et al., 2011), PKP2 [16] (Franchi et al., 2016) and STK38L, ZFHX3, ZNRF1 [17] (Hatzirodos et al., 2014). In sheep, LIG1 and SPP2 [18] (Smith et al., 2019) were differentially expressed between fetal ovaries of fetuses whose dams were exposed to either maintenance or restricted nutrition. GRIK3, EYA2 and BSPH1 were part of a group of genes that were differentially expressed between a subset of Finnsheep and F1 crossbred ewes maintained on a flushing diet [19] (Pokharel et al., 2018). EYA2 was differentially expressed between uniparous and multiparous goat ovaries [20] (Ling et al., 2015). MAP7D1 and TRAPPC3 were associated with NLB and NLW in the Cradock Merino population, while SPP2 was also associated with NLB, supporting the above findings. In the Grootfontein Merino population, the ZFHX3 gene was linked to TWW. The gene SIX6 has been noted as a regulator of gonadotropin releasing hormone (GnRH) in cattle [21] (Cánovas et al., 2014) and sheep [22] (Mellon et al., 2018). In the Afrino population, SIX6 was linked to all three reproductive traits, NLB, NLW and TWW. According to [23] Grive et al. (2014) TAF4B is expressed in both somatic and germ cells in the ovary, and possibly play multiple roles in primordial follicle formation. It is also required for the initial establishment of the primordial follicle reserve at birth. Furthermore, infertility in TAF4B null female mice was associated with defects in early follicle formation and oocyte maturation [24] (Falender et al., 2005a). One of the gene ontology biological processes of TAF4B is “Positive regulation of stem cell proliferation” [25] (https://david.ncifcrf.gov/), which concurs with the associations found in literature. This gene was associated with TWW in the Afrino population in the current study.

The current study only identified one gene that was associated with conception rate in previous literature. PKP2 was one of the genes identified in a genome-wide association study that modulate conception rates in cattle [26] (Sugimoto et al., 2013). A range of genes have been identified that is associated with embryonic development in cattle (PKP2 and IER5L – [27] Killeen et al., 2014), sheep (EYA2 – [28] Ahbara et al., 2019) and mice (CUL7 – [29] Skaar et al., 2005; EYA2 – [30] Grifone et al., 2007; PTK7 – [31] Yen et al., 2009). In the current study, IER5L was linked to NLB in the Cradock Merino population, while PTK7 and CUL7 was associated with TWW in the Afrino population.

Several of the identified genes in the current study were previously associated with the uterus and placenta, as well as maintenance of pregnancy. CABP5 (sheep – [32] Burns et al., 2018), CAP1 (cattle – [33] Fortes et al., 2018), RBP1 and RPB2 (pig – [34] Ma et al., 2018) were associated with the uterus, while CUL7 was involved in mice placental development [35] (Tsunematsu et al., 2006). According to [25] DAVID (https://david.ncifcrf.gov/), one of the GO biological processes of CUL7 is “Placental development”. CABP5 and ELSPBP1 were part of a list of genes in genomic windows that explained more than 1% of the additive genetic variance for early pregnancy in Nelore heifers [36] (Oliveira Júnior et al., 2017). The IL1RN gene was involved in maternal recognition of pregnancy in cattle [37] (Mamo et al., 2012) and early pregnancy in sheep [38] (Pokharel et al., 2020). PKP2 was upregulated in the endometrium during the preimplantation stage of pregnancy in Finnsheep [38] (Pokharel et al., 2020) and, together with STK38L, was expressed in the endometrium of pigs at day 12 of pregnancy [39] (Kim et al., 2012).

Genes previously associated with milk production in cattle were ITFG1 [40] (Mei et al., 2018), ZNF496 [41] (Golik et al., 2011) and TRPS1 [42] (Do et al., 2017), while ZNRF1 was identified as a candidate gene for milk production in Valle del Belice dairy sheep [43] (Sutera, 2018).

Some of the identified genes in this study were also previously linked to spermatogenesis and male reproductive efficiency. TAF4B is important in spermatogonial stem cell maintenance, with high expression in spermatogonia, and some expression in spermatids [44, 45] (Falender et al., 2005b; Cooke et al., 2006). In mice lacking C14orf39, a complete arrest of spermatogenesis, as well as reduced testis size were observed [46] (Gómez-H et al., 2016). In the Afrino population, C14orf39 was associated with all three reproductive traits, NLB, NLW and TWW. RBMS3 was identified as one of the putative genes related to sperm fertility in Assaf sheep rams [47] (Serrano et al., 2019), while ELSPBP1 was differentially expressed in rams with high and low sperm motility [48] (Zhu et al., 2020). ELSPBP1 was associated with dead spermatozoa in cattle [49] (D’Amours et al., 2012) and consequently associated in various studies with spermatozoa in sub-fertile bulls [50, 51] (D’Amours et al., 2010; Kumar et al., 2016). The gene BSPH1 was indicated as a marker of sperm fertility in mice [52] (Heidari-Vala et al., 2020) and linked to sperm capacitation [53, 54] (Plante & Manjunath, 2015; Vala et al., 2018). Apart from being classified under the GO “Development processes” in the Afrino gene DAVID analysis, “Sperm capacitation” is another GO biological process of the BSPH1 gene [25] (https://david.ncifcrf.gov/). The HDAC9 gene was downregulated in 90-day old male goat kids and this was conducive to the vigorous development of spermatogenesis in this period [55] (Bo et al., 2020). TRPM8 is a testosterone receptor [56, 57] (Asuthkar et al., 2015a; Asuthkar et al., 2015b) and TRPM8 channels may be implicated in various physiological processes regulated by androgens [58] (Sutton et al., 2018). TRPM8 was also part of the genes classified under the molecular function “Channel activity” in the Cradock Merino DAVID gene analysis [25] (https://david.ncifcrf.gov/).

Several genes identified in the current study were previously linked with body weight and various growth related traits in sheep, cattle and pigs. The HDAC9 and EYA2 genes were reported to be involved with myogenesis in these species. HDAC9 was associated with myogenesis and muscle development in sheep [59] (Cheng et al., 2020), cattle [60] (De Vos, 2018) and pigs [61] (Zhang et al., 2014). The EYA2 gene was linked to muscle differentiation and development in cattle [62, 63] (Heanue et al., 1999; Hudson et al., 2013), goats [64] (Ling et al., 2019) and pigs [65] (Pérez-Montarelo et al., 2012). EYA2 was also associated with BW in the Grootfontein Merino population. Downregulated expression of the HDAC9 gene has been observed in callipyge animals relative to non-callipyge sheep [66] (Vuocolo et al., 2007). Furthermore, EYA2 was reported to be associated with muscle hypertrophy in goats [64] (Ling et al., 2019).

Some genes were specifically associated with muscle or fat in the live animal. For example, MAP7D1 (cattle – [67] Sweeney et al., 2016), SPP2 (sheep – [68, 69] Trukhachev et al., 2016a; Trukhachev et al., 2016b) and LIG1 (pigs – [70] Da Costa et al., 2004) were associated with or expressed in muscle, while TRAPPC3 (sheep – [71] González-Calvo et al., 2017), THRAP3 (cattle – [72] Perez et al., 2010), PIK3C2A (pigs – [73] Kim et al., 2015) and STK38L (cattle – [74] Lim et al., 2013) were associated with intramuscular or subcutaneous fat. From the genes mentioned here, LIG1 and STK38L were associated with BW in the Afrino population, while XIRP2 was linked to BW and TWW in the Grootfontein Merino population.

Genes identified in the literature to be generally associated with growth are ZFHX3 (cattle – [75] Xu et al., 2017; goats – [76, 77] Zhang et al., 2015; Wei et al., 2018), SIX6 (cattle – [78] Huai et al., 2011; goats – [79] Pan et al., 2011), EYA2 (cattle – [80] Somavilla, 2015) and RBMS3 (cattle – [81] Widmann et al., 2013). Genes associated with body weight at specific ages were the TRPS1 gene which was linked to fetal [82] (Xu et al., 2014) and post weaning weight in sheep [11] (Zhang et al., 2013) and the HDAC9 and LIG1 genes which were associated with birth weight in pigs [61] (Zhang et al., 2014) and cattle [83] (Cole et al., 2014) respectively. Most of these genes were associated with BW in the Afrino or Grootfontein Merino populations. SPP2 was differentially expressed in the loin muscle of Merino sheep with high and low body weight [68, 69] (Trukhachev et al., 2016a; Trukhachev et al., 2016b), while PIK3C2A was associated with body weight in pigs [84] (Bovo et al., 2020). IL1RN was associated with carcass weight in cattle [85] (Daetwyler et al., 2012). Two of the GO biological processes of IL1RN are “Lipid metabolic process” and “Response to glucocorticoid” [25]. This gene was categorised under the “Response to lipid” biological process in the Afrino DAVID gene analysis.

From the results of this study, as well as evidence from other studies performed on sheep and cattle available in literature, some suggestive SNPs and genes with pleiotropic effects were identified that warrant further investigation. The following SNPs were either linked to more than one trait or gene, or the genes associated with these SNPs were previously associated with both reproductive and body weight traits in sheep and cattle. Such a SNP in the Cradock Merino population is OAR1_10554666.1, which is associated with both NLB and NLW and the genes MAP7D1, TRAPPC3 and THRAP3. MAP7D1 is expressed in bovine oocytes [13] (Dickinson, 2016) and bovine muscle [67] (Sweeney et al., 2016), TRAPPC3 in bovine oocytes [14] (Donnison & Pfeffer, 2004) and subcutaneous fat of lambs [71] (González-Calvo et al., 2017), while THRAP3 is involved in the fat profile of bovine muscle [72] (Perez et al., 2010). SNP s20120.1 on OAR1 is associated with NLB and the genes TRMP8, a testosterone receptor [58] (Sutton et al., 2018), and SPP2, previously linked with body weight and expressed in ovine muscle [68, 69] (Trukhachev et al., 2016a; Trukhachev et al., 2016b) and sheep fetal ovaries [18] (Smith et al., 2019). The HDAC9 gene (SNP OAR4_28811142.1), associated with NLB, NLW and TWW in the Cradock Merino population, was previously linked to spermatogenesis in goats [55] (Bo et al., 2020). This gene is involved in myogenesis and muscle development in sheep, cattle and pigs, was observed in callipyge sheep [66] (Vuocolo et al., 2007) and was associated with birth weight in pigs [61] (Zhang et al., 2014).

ZFHX3 linked to OAR14_39202046.1 and associated with TWW, was expressed in bovine follicles [17] (Hatzirodos et al., 2014) and previously linked to growth in cattle [75] (Xu et al., 2017) and goats [77] (Wei et al., 2018). SNP OAR7_76295917.1 in the Afrino population was associated with NLB, NLW and TWW, as well as with two genes, SIX6 and C14orf39. SIX6 is a regulator of GnRH in cattle and sheep [21, 22] (Cánovas et al., 2014; Mellon et al., 2018) and involved in puberty [86] (Fortes et al., 2016) and growth [78] (Huai et al., 2011) in cattle. C14orf39 was previously linked to spermatogenesis [46] (Gómez-H et al., 2016). The TAF4B (SNP OAR23_32551191.1) gene warrants further investigation due to its association with TWW in the Afrino population and evidence from literature linking it to primordial follicle formation [23] (Grive et al., 2014) and sperm stem cell maintenance [45] (Cooke et al., 2006).

SNPs associated with body weight in the current study are OAR9_64654880.1, OAR13_80729511.1 and OAR19_4811675.1 in the Grootfontein Merino population and OAR3_201667351.1 and OAR14_56900862.1 in the Afrino population. The gene TRSP1 (SNP OAR9_64654880.1) was previously associated with milk production in cattle [42] (Do et al., 2017) and fetal and post weaning weight in sheep [11, 82] (Zhang et al., 2013; Xu et al., 2014). SNP OAR13_80729511.1 might be an important marker, due to its associated gene, EYA2 being involved with various reproductive and growth traits. EYA2 is expressed in the ovaries of sheep [19] (Pokharel et al., 2018) and goats [20] (Ling et al., 2015) and involved in embryonic development in sheep [28] (Ahbara et al., 2019). It is also linked to myogenesis [62] (Heanue et al., 1999) and growth [80] (Somavilla, 2015) in cattle. Another gene associated with body weight in the Grootfontein Merino population, RBMS3 (SNP OAR19_4811675.1), was also previously linked to growth in cattle [81] (Widmann et al., 2013) and to sperm fertility in sheep [47] (Serrano et al., 2019). STL38L (SNP OAR3_201667351.1) in the Afrino population was associated with marbling score in cattle [74] (Lim et al., 2013) and expressed in bovine follicles [17] (Hatzirodos et al., 2014) and the porcine endometrium [39] (Kim et al., 2012). The most important SNP associated with body weight in the Afrino population was OAR14_56900862.1, linked to four genes (BSPH1, LIG1, CABP5, ELSPBP1). These genes were associated with sperm fertility in sheep and cattle and pregnancy in sheep and cattle, expressed in sheep ovaries and muscle and involved in birth weight of pigs.

Conclusions

In the current GWAS, 16 SNP markers associated with the reproductive traits were identified among the three populations, while 15 SNPs were associated with body weight. These SNPs were respectively linked to 26 and 21 documented genes in the sheep genome. These genes were previously reported to be associated with reproductive related, as well as with growth related traits in different farm animal species. Continuous long-term selection for both body weight and reproductive traits in these flocks have thus favoured pleiotropic genes influencing both these traits, with the result that SNPs associated with such genes were detected in the association analyses for both groups of traits. Selection for the composite trait TWW most probably especially favoured genes with pleiotropic effects on growth and reproduction.

From the results of this study, supported by other studies performed on sheep and cattle, some genes were identified that warrant further investigation as to their functions and processes relating to growth and reproduction in sheep. These include the genes MAP7D1, TRAPPC3, THRAP3, TRMP8, SPP2, HDAC9, ZFHX3, SIX6, C14orf39, TAF4B, TRSP1 EYA2, RBMS3, STL38L, BSPH1, LIG1, CABP5 and ELSPBP1. Further investigation into these results are necessary, maybe with the inclusion of flocks with different long-term selection objectives and criteria regarding reproduction and body weight, to substantiate these findings. Results will elucidate the metabolic pathways underlying economically important traits, and shed light on the impact of current selection strategies on a genetic level.

Methods

Animals and data

Phenotypic data, pedigrees and genotypes available on animals of the Afrino flock at the Carnarvon Experimental Station, the fine wool Merino stud at Cradock Experimental Station and the Merino stud at Grootfontein Agricultural Development Institute (GADI) were used in this study. These resources were obtained from the GADI-biobank.

The Carnarvon Afrino flock is kept under natural veld conditions at the Carnarvon Departmental Experimental Station (30º 59' S, 22º 9' E) near Carnarvon in the north-western Karoo region of South Africa. The Grootfontein Merino stud is kept at Grootfontein Agricultural Development Institute near Middelburg (31° 28' S, 25° 1' E) in the Eastern Cape Province, while the Cradock fine wool Merino stud is run on irrigated pastures at the Cradock Experimental Station near Cradock (32° 13' S, 28° 41' E) in the Eastern Cape Province.

Selection practices followed in the three flocks since their establishment were discussed in detail by [87] Sűllwald (2020). The Afrino flock has been selected for reproductive traits and body weight since breed development started in 1969 at the Carnarvon Experimental Station. Since 1985, ewes were selected on total weight of lamb weaned, and body weight was one of the main selection objectives of ram selection.

The Grootfontein Merino stud was established at the Grootfontein Agricultural Development Institute in 1955. During the early years, selection was based on subjectively assessed traits. From 1985 onwards the main selection objectives were amended to increase body weight, maintain clean fleece weight and decrease mean fibre diameter and pleat score, while the ewes were selected on total weight of lamb weaned. In 2004, selection objectives were changed, and breeding sires and dams were selected on the basis of the national selection index for relative economical value including reproduction [1] (Herselman & Olivier, 2010).

The Cradock Merino flock was established in 1988 as a genetic fine wool Merino stud at the Cradock Experimental Station. Initially, selection was aimed at improving body weight and maintaining fibre diameter. Until 1999, no selection for reproduction was carried out, while ewes with poor reproductive performance were culled since 2000. As for the Grootfontein Merino flock, selection objectives were changed in 2004 and selection was done on the basis of the national selection index for relative economical value including reproduction [1].

Phenotypic traits included in the study were body weight recorded at selection age at 14 months of age (BW), lifetime number of lambs born (NLB), lifetime number of lambs weaned (NLW) and lifetime total weight of lamb weaned (TWW). Estimated breeding values (EBV) of all traits for the individual animals were obtained as back solutions with the ASReml program [88] (Gilmour et al., 2014). Animal models including direct and maternal additive genetic random effects were fitted for body weight, while only a random direct genetic effect was fitted for the reproductive traits.

For each flock, animals with a range of high and low EBVs for all traits were selected amongst the animals with available genotypes in the GADI-Biobank. Genotypic data were obtained with the Illumina^® Ovine SNP50 BeadChip (Illumina Inc., San Diego, CA). A total of 411 genotypes, comprising 152 Afrino, 129 Cradock Merino and 130 Grootfontein Merino animals, were included.

Genomic analyses

Quality control

The individual genomic datasets for each population were updated with Oar v4.0 SNP Chimp that was downloaded from the SNPchiMp v.3 database [89] (Nicolazzi et al., 2015). Only the 26 autosomal chromosome pairs were used for downstream data analysis. Each dataset’s information was updated for individual identification number, breed, parentage and sex using PLINK v1.07 software [90] (Purcell, 2017). Individual and marker-based quality control (QC) measures were performed on each individual dataset using PLINK v1.07 software [90] (Purcell, 2017). All non-informative SNPs and individuals with missing genotypes were removed at the following parameters: individual call rate of below 90%, a SNP call rate lower than 95%, minor allele frequency of less than 0.02 (MAF <2%) and violation of Hardy-Weinberg equilibrium (p<0.001). Five animals, three from the Grootfontein Merino dataset, and one each from the Afrino and Cradock Merino datasets, were excluded from downstream analyses. The total number of SNPs retained was 42117 for the Afrino, 46196 for the Cradock Merino and 43655 for the Grootfontein Merino datasets.

Principal component analysis (PCA)

Principal Component Analysis was performed to investigate the genetic relatedness of individuals within and between the populations using the merged dataset comprising all three populations. The genomic relationship matrix and estimated principal components were generated with the use of the Genome-wide complex trait analysis v1.24 (GCTA) software [91] (Yang et al., 2011).

Genome-wide association study (GWAS)

The Afrino, Cradock Merino and Grootfontein Merino datasets were analysed separately for each population and each trait, using the software, efficient mixed model association eXpedited (EMMAX) [92] (Kang et al., 2010). EMMAX software was favoured for the analysis as it controls genome-wide error rate successfully compared to other genomic software and it is more suited to calculate the kinship matrix for smaller populations [93] (Eu-ahsunthornwattana et al., 2014).

Results from the association analyses were visualised by creating manhattan plots in R-studio [94] (R Core Team, 2017). Quantile-quantile (Q-Q) plots for each GWAS were also created in R-studio to indicate any deviation from the distribution under the null hypothesis of no association. From the results of the association analyses significant (p<10^-7) and suggestive (p<10^-4) SNP markers were identified on a chromosome-wide level and investigated further on a molecular level. Each SNP marker’s reference SNP name and position were identified from the database Oar v4.0 SNP Chimp that was downloaded from the SNPchiMp v.3 database [89] (Nicolazzi et al., 2015). SNP effects were determined by including all suggestive SNPs for each trait per population in a stepwise regression analysis.

Gene ontology

The Ensembl database [95,96] (Zerbino et al., 2016; Hunt et al., 2018) was used to identify genes that are linked to or in close proximity to the respective SNP markers. To determine the functions of and possible relationships between the identified genes, the genes for each population were analysed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) platform [97,25] (Dennis et al., 2003; https://david.ncifcrf.gov/), with the sheep genome OAR_v4, to categorise the genes in the three gene ontology (GO) categories: biological process, molecular function and cellular component.

Declarations

Ethics approval

Samples for genotyping were collected under approval numbers AP10/3/3 and AP10/3/4 of the Animal Research Ethics Committee of the Grootfontein Agricultural Development Institute. Approval for the use of external data was granted by the ethics committee of the Faculty of Natural and Agricultural Sciences, University of Pretoria (NAS125/2019). I hereby confirm the study was carried out in compliance with the ARRIVE guidelines.

Consent for publication

Not applicable.

Availability of data and materials

The datasets analysed during the current study are not publicly available due it being part of the Grootfontein Biobank resources. The datasets used are, however, available from the corresponding author on reasonable request for specific projects.

Competing interests

The authors declare that they have no competing interests.

Funding

Funding for some of the genotypes was provided by Cape Wools South Africa under project AP1/17/4.

Authors' contributions

MAS - Concept of the study, assist with writing of project proposal, assist with interpretation of results, co-supervisor of the study, writing of the manuscript

SS - Writing of project proposal, analysis of data, interpretation of data for purpose of MSc thesis, assist with writing of the manuscript

WJO - Assist with writing of the manuscript

CV - Concept of the study, assist with writing of project proposal, assist with interpretation of results, supervisor of the study, assist with writing of the manuscript

All the authors have read and approved the final manuscript.

Acknowledgements

The GADI-Biobank, the Eastern Cape Department of Rural Development and Agrarian Reform and the Northern Cape Department of Agriculture, Land Reform and Rural Development are acknowledged for the use of the resources.

References

1. Herselman MJ, Olivier WJ. Description of a model for the calculation of breeding values for profitability. Grootfontein Agric. 2010;10:67–75.
2. Snyman MA, Erasmus GJ, Van Wyk JB, Olivier JJ. Genetic and phenotypic correlations among production and reproduction traits in Afrino sheep. S Afr J Anim Sci. 1998;28(2):74–81.
3. Safari E, Fogarty NM, Gilmour AR. A review of genetic parameter estimates for wool, growth, meat and reproduction traits in sheep. Livest Prod Sci. 2005;92:271–289. doi.org/10.1016/j.livprodsci.2004.09.003.
4. Olivier WJ. The evaluation of a South African fine wool genetic resource flock. PhD thesis, University of Stellenbosch, Stellenbosch, South Africa. 2014. https://doi.org/10.13140/RG.2.1.3436.6960.
5. Zishiri OT, Cloete SWP, Olivier JJ, Dzama K. Genetic parameters for growth, reproduction and fitness traits in the South African Dorper sheep breed. Small Rumin Res. 2013;112:39–48. doi.org/10.1016/j.smallrumres.2013.01.004.
6. Matebesi-Ranthimo PAM, Cloete SWP, Van Wyk JB, Olivier JJ. Genetic parameters for ewe reproduction with objectively measured wool traits in Elsenburg Merino flock. S Afr J Anim Sci. 2017;47:712–721. DOI: 10.4314/sajas.v47i5.15.
7. Xu S-S, Li M-H. Recent advances in understanding genetic variants associated with economically important traits in sheep (Ovis aries) revealed by high-throughput screening technologies. Front Agric Sci Eng. 2017;4:279–288. doi.org/10.15302/JFASE-2017151.
8. Gebreselassie G, Berihulay H, Jiang L, Ma Y. Review on Genomic Regions and Candidate Genes Associated with Economically Important Production and Reproduction Traits in Sheep (Ovis aries). Animals. 2020;10:33. doi:10.3390/ani10010033.
9. Al-Mamun HA, Kwan P, Clark SA, Ferdosi MH, Tellam R, Gondro, C. Genome-wide association study of body weight in Australian Merino sheep reveals an orthologous region on OAR6 to human and bovine genomic regions affecting height and weight. Genet Sel Evol. 2015;47:66. doi.org/10.1186/s12711-015-0142-4.
10. Matika O, Riggio V, Anselme-Moizan M, Law AS, Pong-Wong R, Archibald AL, et al. Genome-wide association reveals QTL for growth, bone and in vivo carcass traits as assessed by computed tomography in Scottish Blackface lambs. Genet Sel Evol. 2016;48:11. doi.org/10.1186/s12711-016-0191-3.
11. Zhang L, Liu J, Zhao F, Ren H, Xu L, Lu J, et al. Genome-wide association studies for growth and meat production traits in sheep. PLoS One. 2013;8(6):e66569. doi.org/10.1371/journal.pone.0066569.
12. Almasi MM, Zamani PP, Mirhoseini SZZ, Moradi MH. Genome wide association study of weaning traits in Lori Bakhtiari sheep. Annals Anim Sci. 2020;20(3):811–824. doi: 10.2478/aoas 2020 0014.
13. Dickinson SE. Effect of pre-ovulatory follicle size on oocyte transcript abundance in beef cows. MSc thesis. University of Missouri. 2016. https://mospace.umsystem. edu/xmlui/bitstream/handle/10355/56052/research.pdf?sequence = 2.
14. Donnison M, Pfeffer PL. Isolation of Genes Associated with Developmentally Competent Bovine Oocytes and Quantitation of Their Levels During Development. Biol Reprod. 2004;71:1813–1821. doi: 10.1095/biolreprod.104.032367.
15. Mamo S, Carter F, Lonergan P, Leal CLV, Naib AA, McGettigan P, et al. Sequential analysis of global gene expression profiles in immature and in vitro matured bovine oocytes: potential molecular markers of oocyte maturation. BMC Genomics. 2011;12: 151. http://www.biomedcentral.com/1471-2164/12/151.
16. Franchi FF, Fontes PK, Santos PH, Razza EM, Loureiro B, Silveira JC, et al. Exosomes of follicular fluid from superstimulated Nelore cows modulate gene expression in bovine cumulus cells and blastocysts. Abstracts A159. In Proceedings 30th Annual Meeting of the Brazilian Embryo Technology Society (SBTE): 26–27 August 2016; Brazil, and 32nd Meeting of the European Embryo Transfer Association (AETE): 9–10 September 2016; Barcelona Spain.
17. Hatzirodos N, Hummitzsch K, Irving-Rodgers HF, Harland ML, Morris SE, Rodgers RJ. Transcriptome profiling of granulosa cells from bovine ovarian follicles during atresia. BMC Genomics. 2014;15:40. www.biomedcentral.com/1471-2164/15/40.
18. Smith P, Juengel J, Maclean P, Rand C, Stanton JL. Gestational nutrition 2: gene expression in sheep fetal ovaries exposed to gestational under nutrition. Reprod. 2019;157(1):13–25. doi.org/10.1530/REP-18-0358.
19. Pokharel K, Peippo J, Honkatukia M, Seppälä A, Rautiainen J, Ghanem N, et al. Integrated ovarian mRNA and miRNA transcriptome profiling characterizes the genetic basis of prolificacy traits in sheep (Ovis aries). BMC Genomics. 2018;19:104. doi: 10.1186/s12864-017-4400-4.
20. Ling YH, Quan Q, Xiang H, Zhu L, Chu MX, Zhang XR, et al. Expression profiles of differentially expressed genes affecting fecundity in goat ovarian tissues. Genet Mol Res. 2015;14(4):18743–18752. dx.doi.org/10.4238/2015.December.28.23.
21. Cánovas A, Reverter A, Deatley KL, Ashley RL, Colgrave ML, Silver GA, et al. Multi-Tissue Omics Analyses Reveal Molecular Regulatory Networks for Puberty in Composite Beef Cattle. PLoS One. 2014;9:1–18. doi.org/10.1371/journal.pone.0102551.
22. Mellon PL, Pandolfi EC, Hoffmann HM. The roles of homeodomain proteins in GnRH neuronal development and fertility. Plenary Presentations - Monday, 16 July 2018. PL05.001 In: Ed H. Dardente. New mechanisms underpinning seasonality in sheep. 9. International Congress of Neuroendocrinology. ICN2018, Jul 2018, Toronto, Canada.
23. Grive KJ, Seymour KA, Mehta R, Freiman RN. TAF4b promotes mouse primordial follicle assembly and oocyte survival. Dev Biol. 2014;392(1):42–51. doi.org/10.1016/j.ydbio.2014.05.001.
24. Falender AE, Shimada M, Lo YK, Richards JS. TAF4b, a TBP associated factor, is required for oocyte development and function. Dev Biol. 2005;288:405–419. doi:10.1016/j.ydbio.2005.09.038.
25. DAVID Bioinformatics Resources 6.8 Database [https://david.ncifcrf.gov/] Accessed on 8 October 2020.
26. Sugimoto M, Sasaki S, Gotoh Y, Nakamura Y, Aoyagi Y, Kawahara T. Genetic variants related to gap junctions and hormone secretion influence conception rates in cows. Proc Nas Academy Sci. 2013;110(48):19495–19500. doi.org/10.1073/pnas.1309307110.
27. Killeen AP, Morris DG, Kenny DA, Mullen MP, Diskin MG & Waters SM. Global gene expression in endometrium of high and low fertility heifers during the mid-luteal phase of the estrous cycle. BMC Genomics. 2014;15:234. doi:10.1186/1471-2164-15-234.
28. Ahbara A, Bahbahani H, Almathen F, Al Abri M, Agoub MO, Abeba A, et al. Genome-Wide Variation, Candidate Regions and Genes Associated With Fat Deposition and Tail Morphology in Ethiopian Indigenous Sheep. Front Genet. 2019;9:699. doi: 10.3389/fgene.2018.00699.
29. Skaar JR, Arai T, DeCaprio JA. Dimerization of CUL7 and PARC Is Not Required for All CUL7 Functions and Mouse Development. Mol Cell Biol. 2005;25(13):5579–5589. doi:10.1128/MCB.25.13.5579–5589.2005.
30. Grifone R, Demignon J, Giordani J, Niro C, Souil E, Bertin F, et al. Eya1 and Eya2 proteins are required for hypaxial somitic myogenesis in the mouse embryo. Dev Biol. 2007;302:602–616. doi:10.1016/j.ydbio.2006.08.059.
31. Yen WW, Williams M, Periasamy A, Conaway M, Burdsal C, Keller R, et al. PTK7 is essential for polarized cell motility and convergent extension during mouse gastrulation. Development. 2009;136:2039–2048. doi:10.1242/dev.030601.
32. Burns GW, Brooks KE, O’Neil EV, Hagen DE, Behura SK, Spencer TE. Progesterone effects on extracellular vesicles in the sheep uterus. Biol Reprod. 2018;98(5):612–622. https://doi.org/10.1093/biolre/ioy011.
33. Fortes MRS, Zacchi LF, Nguyen LT, Raidan F, Weller MMDCA, Choo JJY, et al. Pre- and post-puberty expression of genes and proteins in the uterus of Bos indicus heifers: the luteal phase effect post-puberty. Anim Gen. 2018; 49:539–549. doi: 10.1111/age.12721.
34. Ma X, Li PH, Zhu MX, He LC, Sui SP, Gao S, et al. Genome-wide association analysis reveals genomic regions on Chromosome 13 affecting litter size and candidate genes for uterine horn length in Erhualian pigs. Animal. 2018;12:2453–2461. doi:10.1017/S1751731118000332.
35. Tsunematsu R, Nishiyama M, Kotoshiba S, Saiga T, Kamura T, Nakayama KI. Fbxw8 Is Essential for Cul1-Cul7 Complex Formation and for Placental Development. Mol Cell Biol. 2006;26(16):6157–6169. doi:10.1128/MCB.00595-06.
36. Oliveira Júnior GA, Perez BC, Cole JB, Santana MHA, Silveira J, Mazzoni G, et al. Genomic study and Medical Subject Headings enrichment analysis of early pregnancy rate and antral follicle numbers in Nelore heifers. J Anim Sci. 2017;6(11):4796–4812. https://doi.org/10.2527/jas2017.1752.
37. Mamo S, Mehta JP, Forde N, McGettigan P, Lonergan, P, et al. Conceptus-Endometrium Crosstalk During Maternal Recognition of Pregnancy in Cattle. Biol Reprod. 2012;87(1):6, 1–9. DOI 10.1095/biolreprod.112.099945.
38. Pokharel K, Peippo J, Weldenegodguad M, Honkatukia M, Li M-H, Kantanen J. Gene Expression Profiling of Corpus luteum Reveals Important Insights about Early Pregnancy in Domestic Sheep. Genes. 2020;11(4):415. doi.org/10.3390/genes11040415.
39. Kim M, Seo H, Choi Y, Shim J, Kim H, Lee C-K. Microarray Analysis of Gene Expression in the Uterine Endometrium during the Implantation Period in Pigs. Asian-Aust J Anim Sci. 2012;25(8):1102–1116. dx.doi.org/10.5713/ajas.2012.12076.
40. Mei C, Wang H, Liao Q, Wang L, Cheng G, Wang H, et al. Genetic Architecture and Selection of Chinese Cattle Revealed by Whole Genome Resequencing. Mol Biol Evol. 2018;35(3):688–699. https://doi.org/10.1093/molbev/msx322.
41. Golik M, Glick G, Reicher S, Shirak A, Ezra E, Zeron Y. Differential expression of ruminant ZNF496 variants: Association with quantitative trait locus affecting bovine milk concentration and fertility. J Dairy Sci. 2011;94:2092–2102. doi: 10.3168/jds.2010-3655.
42. Do DN, Bissonnette N, Lacasse P, Miglior F, Sargolzaei M, Zhao X, et al. Genome-wide association analysis and pathways enrichment for lactation persistency in Canadian Holstein cattle. J Dairy Sci. 2017:100:1955–1970. doi.org/10.3168/jds.2016-11910.
43. Sutera AM. Comparison of Genome Wide Association Studies for milk production traits in Valle del Belice dairy sheep. PhD thesis. Universita Degli Studi Di Palermo. 2018.
44. Falender AE, Freiman RN, Geles KG, Lo KC, Hwang K, Lamb DJ, et al. Maintenance of spermatogenesis requires TAF4b, a gonad-specific subunit of TFIID. Genes Dev. 2005;19:794–803. doi: 10.1101/gad.1290105.
45. Cooke PS, Hess RA, Simon L, Schlesser HN, Carnes K, Tyagi G, et al. The transcription factor Ets-Related Molecule (ERM) is essential for spermatogonial stem cell maintenance and self-renewal. Anim Reprod., 2006;3(2):98–107.
46. Gomez-H L, Felipe NM, Sanchez MM, Davies OR, Ramos I, Gracia IT, et al. C14ORF39/SIX6OS1 is a constituent of the synaptonemal complex and is essential for mouse fertility. Nat Commun. 2016;7:13298 doi: 10.1038/ncomms13298.
47. Serrano MM, Ramón M, Calvo JH, Freire F, Vazquez JM, Jiménez MA, et al. Genome-wide association study of sperm traits in Assaf rams. Book of Abstracts, P39. 37th Congr Int Soc Anim Gen. 7–12 July 2019; Leida Spain.
48. Zhu W, Zhang Y, Ren C, Cheng X, Chen J, Ge Z, et al. Identification of proteomic markers for ram spermatozoa motility using a tandem mass tag (TMT) approach. J Proteomics. 2020;210:103438. doi.org/10.1016/j.jprot.2019.103438.
49. D’Amours O, Frenette G, Bordeleau LJ, Allard N, Leclerc P, Blondin P, et al. Epididymosomes transfer epididymal sperm binding protein 1 (ELSPBP1) to dead spermatozoa during epididymal transit in bovine. Biol Reprod. 2012;18:87–94. doi:10.1095/biolreprod.112.100990.
50. D’Amours O, Frenette G, Fortier M, Leclerc P, Sullivan R. Proteomic comparison of detergent-extracted sperm proteins from bulls with different fertility indexes. Reprod. 2010;139:545–556. doi:10.1530/REP-09-0375.
51. Kumar S, Tinson A, Mulligan BP, Ojha S. Gelatin Binding Proteins in Reproductive Physiology. Indian J Microbiol. 2016;56(4):383–393. DOI 10.1007/s12088-016-0618-0.
52. Heidari-Vala H, Sabouhi-Zarafshan S, Prud’homme B, Alnoman A, Manjunath P. Role of Binder of SPerm homolog 1 (BSPH1) protein in mouse sperm-egg interaction and fertilization. Biochem Biophys Res Comm. 2020;527(2):358–364. doi.org/10.1016/j.bbrc.2020.03.127.
53. Plante G, Manjunath P. Murine binder of sperm protein homolog 1: a new player in HDL-induced capacitation. Reprod. 2015;149:367–376. DOI: 10.1530/REP-14-0559.
54. Vala HH. Binder of Sperm protein interference in sperm-egg interaction. MSc thesis. Université de Montréal. 2018. https://papyrus.bib.umontreal.ca/xmlui/handle/1866/20792.
55. Bo D, Jiang X, Liu G, Xu F, Hu R, Wassie T, et al. Multipathway synergy promotes testicular transition from growth to spermatogenesis in early-puberty goats. BMC Genomics. 2020;21:372. doi.org/10.1186/s12864-020-6767-x.
56. Asuthkar S, Elustondo PA, Demirkhanyan L, Sun X, Baskaran P, Velpula KK, et al. The TRPM8 Protein Is a Testosterone Receptor. I. Biochemical evidence for direct TRPM8-testosterone interactions. J Biol Chem. 2015;290(5):2659–2669. DOI 10.1074/jbc.M114.610824.
57. Asuthkar S, Demirkhanyan L, Sun X, Elustondo PA, Krishnan V, Baskaran P, et al. The TRPM8 Protein Is a Testosterone Receptor. II. Functional evidence for an ionotropic effect of testosterone on TRPM8. J Biol Chem. 2015;290(5):2670–2688. DOI 10.1074/jbc.M114.610873.
58. Sutton CM, Ziegler RL, Austin KJ, Alexander BM. Quantitative comparison of TRPM8 positively stained neurons in the hypothalamus and amygdala of rams categorized behaviorally as low or high sexual performers. Transl Anim Sci. 2018;2018.2:S173-S174. doi: 10.1093/tas/txy021.
59. Cheng J, Zhao H, Chen N, Cao X, Hanif Q, Pi L, et al. Population structure, genetic diversity, and selective signature of Chaka sheep revealed by whole genome sequencing. BMC Genomics. 2020; 21:520. doi.org/10.1186/s12864-020-06925-z.
60. De Vos J. A Genome wide association study of carcass traits based on Real Time Ultrasound in South African Nguni cattle. MScAgric thesis. Universiteit of Pretoria, Pretoria, South Africa. 2018. https://link.springer.com/article/10.1186/s12864-020-06925-z.
61. Zhang L, Zhou X, Michal JJ, Ding B, Li R, Jiang Z. Genome wide screening of candidate genes for improving piglet birth weight using high and low estimated breeding value populations. Int J Biol Sci. 2014;10:236–244. doi: 10.7150/ijbs.7744
62. Heanue TA, Reshef R, Davis RJ, Mardon G, Oliver G, Tomarev S, et al. Synergistic regulation of vertebrate muscle development by Dach2, Eya2, and Six1, homologs of genes required for Drosophila eye formation. Genes Dev. 1999;13:3231–3243.
63. Hudson NJ, Lyons RE, Reverter A, Greenwood PL, Dalrymple BP. Inferring the in vivo cellular program of developing bovine skeletal muscle from expression data. Gene Expr Patt. 2013;13:109–125. doi.org/10.1016/j.gep.2013.02.001.
64. Ling Y, Zheng Q, Sui M, Zhu L, Xu L, Zhang Y, et al. Comprehensive Analysis of LncRNA Reveals the Temporal-Specific Module of Goat Skeletal Muscle Development. Int J Mol Sci. 2019;20(16):3950. doi.org/10.3390/ijms20163950.
65. Pérez-Montarelo D, Hudson NJ, Fernández AI, Ramayo-Caldas Y, Dalrymple BP, Reverter A. Porcine Tissue-Specific Regulatory Networks Derived from Meta-Analysis of the Transcriptome. PLoS ONE. 2012;7(9):e46159. doi:10.1371/journal.pone.0046159.
66. Vuocolo T, Byrne K, White J, McWilliam S, Reverter A, Cockett NE, et al. Identification of a gene network contributing to hypertrophy in Callipyge skeletal muscle. Physiol Genomics. 2007;28:253–272.doi:10.1152/physiolgenomics.00121.2006
67. Sweeney T, Lejeune A, Moloney AP, Monahan FJ, Mc Gettigan P, Downey G, et al. The application of transcriptomic data in the authentication of beef derived from contrasting production systems. BMC Genomics. 2016;17:746. doi: 10.1186/s12864-016-2851-7.
68. Trukhachev V, Skripkin V, Kvochko A, Kulichenko A, Kovalev D, Pisarenko S, et al. Correlation between gene expression profiles in muscle and live weight in Dzhalginsky Merino sheep. Rev Colomb Cienc Pecu. 2016;29:188–198. doi: 10.17533/udea.rccp.v29n3a04.
69. Trukhachev V, Belyaev V, Kvochko A, Kulichenko A, Kovalev D, Pisarenko S, et al. Genes expression profiles in the loin muscle of Manych Merino sheep with different live weight. Bulg J Vet Med. 2016;19(1):19–29. DOI: 10.15547/bjvm.927.
70. Da Costa N, Mcgillivray C, Bai Q, Wood JD, Evans G. Nutrient-Gene Interactions Restriction of Dietary Energy and Protein Induces Molecular Changes in Young Porcine Skeletal Muscles. J Nutr. 2004;134:2191–2199. doi.org/10.1093/jn/134.9.2191.
71. Gonzales-Calvo L, Dervishi E, Joy M, Sarto P, Martin-Hernandez R, Serrano M, et al. Genome-wide expression profiling in muscle and subcutaneous fat of lambs in response to the intake of concentrate supplemented with vitamin E. BMC Genomics. 2017;18:92. doi: 10.1186/s12864-016-3405-8.
72. Perez R, Cañón J, Dunner S. Genes associated with long-chain omega-3 fatty acids in bovine skeletal muscle. J Appl Genet. 2010;51(4):479–487.
73. Kim E-S, Ros-Freixedes R, Pena RN, Baas TJ, Estany J, Rothschild MF. Identification of signatures of selection for intramuscular fat and backfat thickness in two Duroc populations. J Anim Sci. 2015;93(7):3292–3302. doi.org/10.2527/jas.2015-8879.
74. Lim D, Lee S-H, Kim N-H, Cho Y-M, Chai H-H, Seong H-H, et al. Gene Co-expression Analysis to Characterize Genes Related to Marbling Trait in Hanwoo (Korean) Cattle. Asian-Aust J Anim Sci. 2013;26(1):19–29. dx.doi.org/10.5713/ajas.2012.12375.
75. Xu H, Zhang S, Zhang X, Dang R, Lei C, Chen H, et al. Evaluation of novel SNPs and haplotypes within the ATBF1 gene and their effects on economically important production traits in cattle. Arch Anim Breed. 2017;60:285–296. doi.org/10.5194/aab-60-285-2017.
76. Zhang X, Wu X, Jia W, Pan C, Li X, Lei C, et al. Novel Nucleotide Variations, Haplotypes Structure and Associations with Growth Related Traits of Goat AT Motif-Binding Factor (ATBF1) Gene. Asian-Aust J Anim Sci. 2015;28(10):1394–1406. dx.doi.org/10.5713/ ajas.14.0860.
77. Wei Z, Wang K, Hui Y, Yan H, Zhu H, Qu L, et al. Detection of insertion/deletions (indels) of the ATBF1 gene and their effects on growth-related traits in three indigenous goat breeds. Arch Anim Breed. 2018;61;311–319. doi.org/10.5194/aab-61-311-2018.
78. Huai YT, Lan XY, Ma L, Wang J, Lei CZ, Chen H, et al. Novel Mutation in TGA Stop Codon of Bovine SIX6. Gen Mol Biol. 2011;45:190–196. doi.org/10.1134/S0026893310061093.
79. Pan CY, Lan XY, Zhao HY, Hu SR, Huai YT, Lei CZ. A novel genetic variant of the goat SIX6 gene and its association with production traits in Chinese goat breeds. Genet Mol Res. 2011;10:3888–3900.
80. Somaville AL. Prediction of genomic-enabled breeding values and genome-wide association study for feedlot average daily weight gain in Nelore cattle. PhD thesis, Universidade Estadual Paulista. UNESP. 2015. https://repositorio.unesp.br/handle/11449/128158.
81. Widmann P, Reverter A, Fortes MAS, Weikard R, Suhre K, Hammon H, et al. A systems biology approach using metabolomic data reveals genes and pathways interacting to modulate divergent growth in cattle. BMC Genomics. 2013;14:798. http://www.biomedcentral.com/1471-2164/14/798.
82. Xu L, Zhao F, Ren H, Li L, Lu J, Liu J, et al. Co-Expression Analysis of Fetal Weight-Related Genes in Ovine Skeletal Muscle during Mid and Late Fetal Development Stages. Int J Biol Sci. 2014;10(9):1039–1050. doi: 10.7150/ijbs.9737.
83. Cole JB, Waurich B, Bickhart DM, Swalve HH. A genome-wide association study of calf birth weight in Holstein cattle using single nucleotide polymorphisms and phenotypes predicted from auxiliary traits. J Dairy Sci. 2014;97;3156–3172. doi.org/10.3168/jds.2013-7409.
84. Bovo S, Ribani A, Muñoz M, Alves E, Araujo JP, Bozzi R. Whole-genome sequencing of European autochthonous and commercial pig breeds allows the detection of signatures of selection for adaptation of genetic resources to different breeding and production systems. Genet Sel Evol. 2020;52:33. doi.org/10.1186/s12711-020-00553-7
85. Daetwyler HD, Swan AA, Van der Werf JHJ, Hayes BJ. Accuracy of pedigree and genomic predictions of carcass and novel meat quality traits in multi-breed sheep data assessed by cross-validation. Genet Sel Evol. 2012;44:33. http://www.gsejournal.org/content/44/1/33.
86. Fortes MRS, To L, Porto LR, Reverter A, Moore SS, Lehnert SA, et al. Polymorphisms and genes associated with puberty in heifers. Theriogenology. 2016;86:333–339. doi.org/10.1016/j.theriogenology.2016.04.046.
87. Sűllwald S. A genome wide association study of body weight and reproduction traits in two South African sheep breeds. MSc. Thesis. University of Pretoria, Pretoria, South Africa. 2020.
88. Gilmour AR, Gogel BJ, Cullis BR, Welham SJ, Thompson R. ASReml User Guide Release 4.1 Structural Specification, VSN International Ltd, Hemel Hempstead, HP1 1ES, UK. . 2014.
89. Nicolazzi EL, Caprera A, Nazzicari N, Cozzi P, Strozzi F, Lawley C, et al. SNPchiMp v.3: Integrating and standardizing single nucleotide polymorphism data for livestock species. BCM Genomics. 2015;16:283. doi.org/10.1186/s12864-015-1497-1.
90. Purcell S, Neale B, Todd-Brown K. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–575.
91. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: A tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88:76–82. https://doi.org/10.1016/j.ajhg.2010.11.011.
92. Kang HM, Sul JH, Service SK, Zaitlen NA, Kong SY, Freimer NB, et al. Variance component model to account for sample structure in genome-wide association studies. Nat Genet. 2010;42:348–366.
93. Eu-ahsunthornwattana J, Miller EN, Fakiola M, Jeronimo SMB, Blackwell JM, Cordell HJ. Comparison of Methods to Account for Relatedness in Genome-Wide Association Studies with Family-Based Data. PLoS Genet. 2014;10;1–20.
94. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/. 2017.
95. Zerbino DR, Johnson N, Juetteman T, Sheppard D, Wilder SP, Lavidas I, et al. Ensembl regulation resources. Database 2016, 1–13. https://doi.org/10.1093/database/bav119.
96. Hunt SE, McLaren W, Gil L, Thormann A, Schuilenburg H, Sheppard D, et al. Ensembl variation resources. Database Volume 2018. doi:10.1093/database/bay119.
97. Dennis G, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, et al. DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol. 2003;4:R60. https://doi.org/10.1186/gb-2003-4-9-r60.