13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License FULL PAPER doi: 10.1111/age.13113 Detecting the footprint of selection on the genomes of Murciano-Granadina goats D. Guan* , A. Martı́nez†, M. G. Luigi-Sierra* , J. V. Delgado†, V. Landi†,‡, A. Castelló*,§, J. Fernández Álvarez¶, X. Such**, J. Jordana§ and M. Amills*,§ *Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Universitat Autònoma de Barcelona, Bellaterra 08193, Spain. † Departamento de Genética, Universidad de Córdoba, Córdoba 14071, Spain. ‡Department of Veterinary Medicine, University of Bari “Aldo Moro”, SP. 62 per Casamassima km. 3, Valenzano 70010, Italy. §Departament de Ciència Animal i dels Aliments, Facultat de Veterinària, Universitat Autònoma de Barcelona, Bellaterra 08193, Spain. ¶Asociación Nacional de Criadores de Caprino de Raza Murciano-Granadina (CAPRIGRAN), Fuente Vaqueros, Granada 18340, Spain. **Group of Research in Ruminants (G2R), Department of Animal and Food Science, Universitat Autònoma de Barcelona (UAB), Bellaterra 08193, Spain. Summary Artificial selection is one of the major forces modifying the genetic composition of livestock populations. Identifying genes under selection could be useful to elucidate their impact on phenotypic variation. We aimed to identify genomic regions targeted by selection for dairy and pigmentation traits in Murciano-Granadina goats. Performance of a selection scan based on the integrated haplotype score test in a population of 1183 Murciano-Granadina goats resulted in the identification of 77 candidate genomic regions/SNPs. The most significant selective sweeps mapped to chromosomes 1 (69.86 Mb), 4 (41.80–49.95 Mb), 11 (65.74 Mb), 12 (31.24 and 52.51 Mb), 17 (34.76–37.67 Mb), 22 (31.75 Mb), and 26 (26.69–31.05 Mb). By using previously generated RNA-Seq data, we built a catalogue of 6414 genes that are differentially expressed across goat lactation (i.e. 78 days post-partum, early lactation; 216 days post-partum, late lactation; 285 days post-partum, dry period). Interestingly, 183 of these genes mapped to selective sweeps and several of them display functions related with lipid, protein, and carbohydrate metabolism, insulin signaling, cell proliferation, as well as mammary development and involution. Of particular interest are the CSN3 and CSN1S2 genes, which encode two major milk proteins. Additionally, we found three pigmentation genes (GLI3, MC1R, and MITF) co-localizing with selective sweeps. Performance of a genome-wide association study and Sanger sequencing and TaqMan genotyping experiments revealed that the c.801C>G (p.Cys267Trp) polymorphism in the melanocortin 1 receptor (MC1R) gene is the main determinant of the black (GG or GC genotypes) and brown (CC genotypes) colorations of Murciano-Granadina goats. Keywords Casein, coat color, iHS test, MC1R, selective sweeps selection (Harris et al. 2018). Multiple statistical methods Introduction based on diversity reduction, shifts of the site frequency Selection for production and morphological traits has spectrum, specific LD patterns, and single site or haplotype strongly modified the genetic composition of livestock differentiation have been developed in order to detect hard populations (Saravanan et al. 2020). Increased linkage selective sweeps at a genome-wide scale (Vitti et al. 2013; disequilibrium (LD) and reduced haplotype diversity around Saravanan et al. 2020). In livestock, the identification of the selected site are typical genomic signatures left by genes under selection should be helpful to elucidate their biological roles as well as to pinpoint mutations with beneficial effects on traits of economic interest. Aiming to Address for correspondence achieve this goal, multiple selection scans have been M. Amills, Centre for Research in Agricultural Genomics (CRAG), CSIC- performed in goats to map selective sweeps related to IRTA-UAB-UB, Universitat Autònoma de Barcelona, Bellaterra 08193, pigmentation (Bertolini et al. 2018; Oget et al. 2018), milk Spain. traits (Bertolini et al. 2018; Oget et al. 2018; Zhang et al. Email: marcel.amills@uab.cat 2018), litter size (Lai et al. 2016; E et al., 2019), cashmere Accepted for publication 16 June 2021 © 2021 The Authors. Animal Genetics published by 683 John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 684 Guan et al. production (Li et al. 2017; Zhang et al. 2018), and rate (<90%), low frequency (<1%), or that deviated environmental adaptation (Kim et al. 2016; Song et al. significantly from the Hardy–Weinberg expectation 2016; Brito et al. 2017; Onzima et al. 2018), to mention a (P ≤ 1 × 10−5). We also removed samples with genotyping few. rates lower than 90%. After applying these filtering criteria, Murciano-Granadina (MG) goats, which can display 47 098 SNPs and 1183 MG goats were selected to perform black or brown coats, are fundamentally devoted to milk a selection scan. production and cheese manufacturing (Delgado et al. 2017). The breeding program is currently managed by Detecting signatures of selection in the MG breed the Asociación Nacional de Caprino de Raza Murciano- Granadina (CAPRIGRAN) and aims to improve milk yield To detect selective sweeps in the MG population, we and composition as well as morphology (Delgado et al. employed the integrated haplotype score (iHS) algorithm 2017). So far, selection has substantially improved the milk implemented in the REHH 2.0 software (Gautier et al. 2017). productivity of MG goats, which currently reaches 530 kg The iHS algorithm measures the observed decay of extended of milk yield per lactation and 5.6% and 3.6% of fat and haplotype homozygosity summed over both directions away protein contents respectively (Delgado et al. 2017). The from a given SNP (Voight et al. 2006). Large negative iHS main aim of the current work was to identify genomic scores indicate that a derived allele has swept up in signatures produced by artificial selection for dairy and frequency, while large positive iHS values at nearby SNPs pigmentation traits in MG goats. To achieve this goal, we would imply either that the ancestral allele hitchhiked with have genotyped a population of MG goats with the Goat the selected site or that selection switched to favor an SNP50 BeadChip and subsequently, we have investigated ancestral allele (Voight et al. 2006). In this study, we the footprint of selection by integrating several sources of assumed as ancestral the reference allele by using the information. ‘polarize_vcf = FALSE’ option implemented in the REHH package (Gautier et al. 2017). Genome-wide iHS values were adjusted to make sure that they approximately follow Materials and methods a normal distribution with mean 0 and variance 1 (Voight et al. 2006; Gautier & Vitalis 2012; Gautier et al. 2017). Sampling of MG goats and genotyping with the Goat Following Gautier et al. (2017), statistical significance was SNP50 BeadChip inferred as follows: The experimental design of this research is summarized in piHS ¼ log10 ð1  2jΦðiHSÞ  0:5jÞ: Fig. S1. As animal material, we have used 1197 MG goats distributed in 19 farms from the autonomous region of where Φ(iHS) represents the Gaussian cumulative distribu- Andalusia. The current census of the MG breed accounts for tion function and piHS is a two-sided P-value, in a −log10 more than 100 000 individuals, which are generally raised scale, associated with the null hypothesis (H0 = selective in semi-intensive conditions with the goal of producing neutrality). The false discovery rate (FDR) approach was cheese (Spanish Ministry of Agriculture, Fisheries and Food, used to correct for multiple testing, and the significance https://www.mapa.gob.es/es/ganaderia/temas/zootecnia/ threshold was set to a q-value ≤0.05 (Benjamini & razas-ganaderas/razas/catalogo-razas/caprine/murciano- Hochberg 1995). Following Guan et al. (2020), significant granadina/iframe-ejemplo-arca.aspx). One of the main SNPs were grouped into genomic regions when their hallmarks of the MG breed is its extraordinary adaptation physical positions were <1 Mb apart. Genes located within to harsh climate conditions and aridity, being able to thrive regions under selection were retrieved with the BEDTOOLS in marginal lands with scarce food resources (Spanish software (Quinlan & Hall 2010) by using the ARS1 goat Ministry of Agriculture, Fisheries and Food). Blood samples genome (Bickhart et al. 2017) as a reference. Functional were collected in EDTA K3 coated vacuum tubes and stored enrichment analysis for the set of genes co-localizing with at −20°C before processing. Genomic DNA was isolated putative selective sweeps was carried out with the DAVID using a previously reported salting-out procedure (Guan Bioinformatics Resources 6.8 database (Huang et al. 2008, et al. 2020) and resuspended in 1 ml TE buffer (Tris-HCl 2009) by considering goat background genes. 10 mmol/l, EDTA 1 mmol/l, pH = 8). The Goat SNP50 BeadChip (Tosser-Klopp et al. 2014), which contains Establishing a catalogue of genes related with goat 53 347 single nucleotide polymorphisms (SNPs), was used lactation and pigmentation and analyzing their co- to genotype samples following the instructions of the localization with selective sweeps manufacturer (Illumina). Data normalization and genotype calling were performed with the GENOMESTUDIO software 2.0.4 Milk yield and composition are complex traits determined by (Illumina, https://emea.illumina.com). The PLINK v1.9 soft- the action of multiple genes, the identities of which are ware (Purcell et al. 2007) was used to filter out unmapped mostly unknown. To facilitate the interpretation of the and non-autosomal SNPs, as well as those with a low call selection scan data, we have established a catalogue of © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Selective sweeps in Murciano-Granadina goats 685 genes related with lactation. This goal has been achieved by random effects with a n-dimensional multivariate normal reanalyzing RNA sequencing (RNA-Seq) data generated in a distribution (0, λτ−1K), being τ−1 the variance of the previous work (Guan et al. 2020). In the study of Guan et al. residual error, λ the ratio between the two variance (2020), the differential expression of mRNAs transcribed in components, and K a n × n relatedness matrix derived from the mammary gland of goats in the beginning (T1, day 78) 42 793 valid SNPs. Lastly, ϵ represents a vector of errors. In and ending (T2, day 216) of lactation as well as in the dry this model, we did not include the fixed factors defined by period (T3, day 285) was investigated by setting, as Guan et al. (2020) because they are not expected to have thresholds of significance, a q-value ≤0.05 and an absolute any effect on coat color. The GEMMA software contrasts the fold-change >2.82. In the current work, we have used less alternative hypothesis (H1: β ≠ 0) against the null hypoth- stringent thresholds of significance, i.e. q ≤ 0.05 and an esis (H0: β = 0) by carrying out likelihood ratio tests for fold-change >1.5, in order to increase the chances of each marker. Besides, the relatedness matrix, which is finding differentially expressed (DE) genes co-localizing with constructed by taking into account all genome-wide SNPs, selective sweeps (Fig. S1). Briefly, we have filtered out is employed to correct for population structure. Multiple features with read counts below 10 based on the expression testing was implemented through an FDR approach (Ben- matrix of protein-coding genes generated by the FEATURE- jamini & Hochberg 1995), and a q-value ≤0.05 was COUNTS v1.5.0-p1 tool (Liao et al. 2014). Subsequently, established as threshold of significance in the GWAS. differential expression analysis was carried out with the DESEQ2 software (Love et al. 2014). After correcting for Sanger sequencing of the MC1R gene multiple testing with the FDR approach (Benjamini & Hochberg 1995), a list of DE genes was generated by taking Two pairs of primers were designed with the Primer-BLAST into account the thresholds of significance mentioned tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/inde above. Moreover, we have built a catalogue of 80 genes x.cgi) in order to amplify the coding region of the MC1R with well-known roles in pigmentation by consulting the gene (Table S1) in 9 black and 13 brown MG individuals. Color Genes database (http://www.espcr.org/micemut/, PCRs contained 50 ng genomic DNA, 1× BIOTAQ PCR accessed 27 February 2020). Once these two catalogues buffer (Bioline), 200 µmol/l of dNTPs, 0.2 µmol/l of each of genes associated with lactation and pigmentation were primer, 1.5 mmol/l MgCl2, and 0.65 units of BIOTAQ DNA established, we investigated their overlap with the sets of Polymerase (Bioline). Nuclease-free water was added to a genes mapping to selective sweeps (please see the previous final volume of 25 µl. The thermal cycle was as follows: a section). hot-start step at 95°C for 2 min, followed by 34 cycles of 95°C for 45 s (denaturation), 60°C for 45 s (annealing), and 72°C for 45 s (extension), plus a final extension at 72°C Performance of a genome-wide association study for for 5 min. Five µl of the PCR were mixed with 1.5 µl of a pigmentation traits in MG goats mixture containing 1.13 µl PCR buffer 1× (composition for The selection scan highlighted the existence of several 1 ml: 100 µl PCR Gold Buffer 10× + 100 µl MgCl2 pigmentation genes co-localizing with selective sweeps. 25 mmol/l + 800 µl H2O), 0.12 µl Exonuclease I Following Bâlteanu et al. (2021), we sought to confirm (20 units/µl, Thermo Fisher Scientific) and 0.25 µl FastAP whether any of these genes is involved in the coloration of Thermosensitive Alkaline Phosphatase (1 unit/µl, Thermo MG goats by carrying out a genome-wide association study Fisher Scientific). This mixture was incubated at 37°C for (GWAS). Coat color (black or brown) of 529 MG goats 15 min plus 85°C for 15 min. Purified amplicons were chosen at random was phenotyped by visual assessment sequenced with the BIGDYE TERMINATOR CYCLE SEQUENCING KIT v1.1 (Fig. S1). These 529 goats also had Goat SNP50 BeadChip (Applied Biosystems). Sequencing reactions were run on an genotypes. Unmapped and non-autosomal SNPs, as well as ABI 3730 DNA analyzer (Applied Biosystems). Finally, those with a low call rate (<90%) and low frequency (<1%) sequences were viewed and aligned with the MOLECULAR were filtered out with the PLINK v1.9 software (Purcell et al. EVOLUTIONARY GENETICS ANALYSIS software (MEGA X; Kumar 2007). After applying these filtering criteria, 43 240 SNPs et al. 2018). were selected to perform a GWAS analysis with the GEMMA software (version 0.98.1), which implements the Genome- wide Efficient Mixed Model Association algorithm (Zhou & TaqMan genotyping of a missense mutation in the Stephens 2012). The following statistical model was used: MC1R gene To further confirm the causality of the c.801C>G muta- Y ¼ xβ þ u þ ɛ: tion, we used a Custom TaqMan SNP Genotyping Assay where Y is a vector of phenotypic values coded as 1 (black, (Applied Biosystems) to genotype 49 black and 41 brown N = 391) or 2 (brown, N = 138); x is a n-vector of marker individuals. These 90 individuals were randomly selected genotypes harbored by each individual; β is the effect size of from the population used in the GWAS analysis (see the marker (allele substitution effect); u is a n-vector of above). The description of TaqMan probes is reported in © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 686 Guan et al. Table S2. Five samples with GG (N = 1), CC (N = 2), and 17 (34.76–37.67 Mb, iHS = 5.34, q = 6.38 × 10−5). Inter- GC (N = 2) genotypes ascertained by Sanger sequencing estingly, the genomic region on chromosome 6 (86.08– were used as positive controls. Genotyping reactions were 86.21 Mb) overlaps the casein gene cluster (Table S3). The carried out in a final volume of 15 μl containing 1× most significant SNPs on chromosomes 16 and 17 map close TaqMan Universal PCR Master Mix (Applied Biosystems), to the solute carrier family 45 member 1 (SLC45A1) gene and 1× TaqMan Custom Genotyping Assay designed for to the transient receptor potential cation channel subfamily C rs669694251 (Applied Biosystems) and 18 ng genomic member 3 (TRPC3) gene respectively (Table S3). In addition, DNA. Real-time PCRs were performed in 96-well reaction we found selective sweeps on chromosomes 18 (16.64– plates and they were run in a 7900-HT Real Time PCR 19.46 Mb, iHS = −5.76, q = 7.19 × 10−6), and 22 system (Applied Biosystems). The thermal profile was: 50°C (31.75 Mb, iHS = −4.81, q = 7.45 × 10−4), as shown in for 2 min, 95°C for 10 min, and 40 cycles of 95°C for 15 s Fig. 1, Tables 1 & S3. As indicated in Table S4, functional and 60°C for 1 min. Genotypes were obtained using the enrichment analysis of genes mapping to selective sweeps did Genotyping Analysis Module of the APPLIED BIOSYSTEMS not identify significant pathways/terms after correction for ANALYSIS Software accessible in the Thermofisher Cloud multiple testing (q > 0.05), but the pathway ‘oocyte meiosis’ (https://www.thermofisher.com/es/es/home/digital-scie was the most significant one at the nominal level nce.html, accessed March 20, 2020). (P = 8.98 × 10−4). Results Co-localization between selective sweeps and genes differentially expressed across goat lactation Identification of putative selective sweeps in a Selection scans often generate long list of genes mapping to population of MG goats selective sweeps. Such data can be difficult to interpret By making a selection scan based on the iHS score algorithm, because the functional relationship between genes and we identified 77 genomic regions/SNPs that reached the phenotypes under selection is usually unknown. To shed threshold of significance (Fig. 1, Tables 1 & S3). Among some light into this issue, we have built a catalogue of genes them, we observed a highly significant signal on chromosome that are expressed in the caprine mammary gland and DE 4 (41.80–49.95 Mb, iHS = 10.30, q = 3.38 × 10−20), across lactation. By reanalyzing RNA-Seq data generated by which co-localized with 72 genes (Table S3). In this selective Guan et al. (2020), we were able to identify 410 (T1/T2), sweep, the most significant SNP (rs268259933) lies close to 5368 (T1/T3), and 4437 (T2/T3) DE genes. In total, 6414 the zinc finger protein 804B (ZNF804B) gene, and there are genes showed differential expression across comparisons, also several genes related with metabolism, such as adenylyl and 183 of them co-localized with selective sweeps cyclase type 1 (ADCY1), insulin-like growth factor-binding (Table S5). In this reduced list of 183 genes, we found loci proteins 1 (IGFBP1) and 3 (IGFBP3). As shown in Fig. 1 and related with: lipid (abhydrolase domain-containing protein Table 1, three additional selective sweeps of considerable size 2, ABHD2; perilipin 1, PLIN1; stearoyl-CoA desaturase, were observed on chromosome 6 (86.08–86.21 Mb, iHS = SCD), casein (casein αS2, CSN1S2; casein κ, CSN3), and −5.44, q = 4.08 × 10−5), chromosome 16 (40.09– carbohydrate (glucokinase, GCK; enolase 1, ENO1) meta- 51.83 Mb, iHS = 5.23, q = 1.05 × 10−4), and chromosome bolism; insulin signaling (phosphatidylinositol-4-phosphate Figure 1 Selection scan in a population of Murciano-Granadina goats. The x- and y-axes indicate genomic coordinates and negative log10 P-values respectively. The dashed line indicates the threshold of significance (q ≤ 0.05). Significant SNPs are marked in red. © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Selective sweeps in Murciano-Granadina goats 687 Table 1 List of selective sweeps identified with the iHS-based algorithm in Murciano-Granadina goats (the full list of selective sweeps is provided in Table S3). Chromosome Region (Mb) SNP iHS value P q Key genes 4 41.80–49.95 rs268259933 10.3 7.32E-25 3.38E-20 ZNF804B, ADCY1, IGFBP1, IGFBP3, GCK, OGDH, HUS1, GLI3 5 90.12 rs268282714 4.69 2.72E-06 1.19E-03 PIK3C2G 6 86.08–86.21 rs268293092 −5.44 5.47E-08 4.08E-05 CSN3, CSN1S2 7 60.11–60.89 rs268262914 −4.44 8.86E-06 3.39E-03 EGR1, CDC25C 16 40.09–51.83 rs268275948 5.23 1.68E-07 1.05E-04 SLC45A1, ENO1, PIK3CD, PGD, B3GALT6, TNFRSF8, TNFRSF9, TNFRSF18 17 34.76–37.67 rs268275466 5.34 9.24E-08 6.38E-05 TRPC3 18 16.64–19.46 rs268267832 −5.76 8.24E-09 7.19E-06 MC1R 21 5.85–8.98 rs268257442 4.11 3.89E-05 1.12E-02 IGF1R 21 19.88–21.86 rs268283170 4.65 3.27E-06 1.40E-03 ANPEP, ABHD2, PLIN1 22 31.75 rs268246615 −4.81 1.51E-06 7.45E-04 MITF 25 31.54–35.43 rs268242461 −3.91 9.33E-05 2.22E-02 MDH2 26 26.69–31.05 rs268273064 −6.68 2.40E-11 4.11E-08 SCD, NFKB2 3-kinase catalytic subunit type 2 γ, PIK3C2G; part of the P-values observed in the GWAS deviate from the phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic expected P-values drawn from a theoretical Chi-square subunit δ, PIK3CD); Krebs cycle and pentose phosphate distribution (Fig. 2d). pathway (malate dehydrogenase 2, MDH2; oxoglutarate Based on these results, we decided to sequence the coding dehydrogenase, OGDH; phosphogluconate dehydrogenase, region of the MC1R gene. By doing so, we identified two PGD); cell proliferation (cell division cycle 25C, CDC25C; missense mutations: c.748G>T (p.Val250Phe) and HUS1 checkpoint clamp component, HUS1); and mammary c.801C>G (p.Cys267Trp). These two polymorphisms have development and involution (alanyl aminopeptidase, mem- been already registered in the Ensembl database (https:// brane, ANPEP; early growth response 1, EGR1; insulin like www.ensembl.org) with identifiers rs657434682 and growth factor 1 receptor, IGF1R; insulin like growth factor rs669694251 respectively. Visual inspection of the binding protein 3, IGFBP3; nuclear factor kappa B subunit sequences revealed a single G peak (GG genotype) at 2, NFKB2; TNF receptor superfamily members 8, 9 and 18, position 748 in nine black and 10 brown goats, while only TNFRSF8, TNFRSF9 and TNFRSF18). three brown goats showed a double GT peak at this location (Table S6). Moreover, no goat displayed a single T peak (TT genotype). Altogether, these results rule out c.748G>T Genomic architecture of coat color in MG goats (p.Val250Phe) as a causal mutation for the black/brown With regard to coat color, we generated a catalogue of 80 coat color of MG goats. With regard to c.801C>G, black genes with well-known functions on pigmentation (http:// goats showed either a single G peak (GG genotype, N = 1) www.espcr.org/micemut/, accessed 27 February 2020) and or a GC double peak (GC genotype, N = 8), while all brown explored their co-localization with selective sweeps. We goats (N = 13) displayed a single C peak (CC genotype, found that the GLI family zinc finger 3 (GLI3), melanocortin Fig. 3a, Table S6). 1 receptor (MC1R) and melanocyte-inducing transcription By performing a TaqMan genotyping experiment in 49 factor (MITF) pigmentation genes map to selective sweeps black and 41 brown individuals, we found that all brown on chromosomes 4 (41.80–49.95 Mb), 18 (16.64– goats were CC and all black goats were either GG or GC 19.46 Mb) and 22 (31.75 Mb) respectively. To further (Fig. 3b), indicating that the G-allele has probably a clarify the genomic architecture of pigmentation in MG dominant effect over the C-allele. Moreover, a second goats, we performed a GWAS analysis comprising 529 MG GWAS was carried out by using a SNP data set comprising goats with available color phenotypes (Fig. 2a,b). We 43 240 SNPs contained in the Goat SNP50 BeadChip plus detected 27 SNPs on chromosome 18 (12.18–22.30 Mb) the rs669694251 (c.801C>G) marker. The GEMMA software showing significant associations with coloration (Fig. 2c), (Zhou & Stephens 2012) was employed to investigate the while no significant signal was observed for chromosomes 4 association between these 43 241 markers and color and 22. The chromosome 18 marker showing the highest phenotypes in 49 black and 41 brown MG goats (N = 90) statistical significance was rs268287597 (chromosome 18: genotyped with the TaqMan assay. In Fig. 3c, it can be 15 887 136 bp, q = 1.34 × 10−19, Fig. 2c), which lies appreciated that marker rs669694251 displays the most close to the MC1R gene (chromosome 18: 16 104 986– significant association with coat color (c.801C>G, 16 105 939 bp). The quantile-quantile plot illustrates that q = 2.91 × 10−25). In this subset of 90 individuals, we © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 688 Guan et al. Figure 2 Murciano-Granadina goats with black (a) and brown (b) pigmentation patterns, which were respectively coded as 1 and 2 in the genome- wide association study. (c) Manhattan plot depicting the associations between coat color and Goat SNP50 BeadChip genotypes from 391 black and 138 brown Murciano-Granadina goats. Negative log10 P-values (y-axis) of the associations between SNPs and phenotypes are plotted against the genomic location of each SNP marker (x-axis). Markers on different chromosomes are denoted with different colors. The horizontal dashed line indicates statistical significance after correction for multiple testing by using the false discovery rate approach reported by Benjamini & Hochberg (1995). The arrow indicates the leading SNP that shows the highest association with phenotype (rs268287597, q = 1.34 × 10−19). (d) Quantile- quantile plot of the expected versus observed P-values in the genome-wide association study analysis for pigmentation. also found that the rs669694251 SNP is in high disequi- dehydrogenase, an enzyme participating in the Krebs cycle, librium with the leading SNP rs268287597 identified in the one of the main pathways by which acetyl-CoA derived first GWAS analysis (r2 = 0.30, D0 = 0.81). from lipids, sugars, and proteins is oxidized to generate NADH, GTP, and FADH2 as well as the building blocks needed for macromolecule synthesis (Martı́nez-Reyes & Discussion Chandel 2020). On chromosome 4 (41.80–49.95 Mb), we also detected the IGFBP3 gene, which modulates the Selective sweeps encompass genes related with involution of the mammary gland (Neuenschwander et al. mammary metabolism and lactation 1996; Flint et al. 2001) as well as lipid and carbohydrate The most significant selective sweep (q = 3.38 × 10−20) metabolism (Kim 2013). mapped to chromosome 4 (41.80–49.95 Mb) and harbored The iHS scan also captured a significant signal 72 genes. The leading SNP was located downstream of the (q = 4.08 × 10−5) on a chromosome 6 region (86.08– ZNF804B gene, which has been associated with dry matter 86.21 Mb) co-localizing with the CSN1S2 and CSN3 genes. intake (de Oliveira et al. 2014) and milk oleic acid content in Casein genes synthesize the major fraction (>80%) of the cattle (Mychaleckyj et al. 2018), as well as with fat milk proteome and their variability has been associated with deposition in sheep (Xu et al. 2017). Moreover, in this milk composition in MG (Guan et al. 2020) and Alpine and chromosome 4 region, we detected several genes that are Saanen (Martin et al. 2017) goats. By using a candidate DE in lactating and non-lactating goats. For instance, the gene approach, Caravaca et al. (2008, 2009, 2011) found GCK gene encodes an enzyme that catalyzes the formation that the BB and EE casein αS1 (CSN1S1) genotypes are of glucose-6-phosphate which happens to be the rate- associated with CSN1S1 protein synthesis and curdling rate limiting step in glucose metabolism (Matschinsky & Wilson respectively. Moreover, the AB and BB CSN3 genotypes 2019). This molecule is a primary glucose sensor in have been associated with higher milk protein and casein mammals, being responsible for the regulation of glucose contents (Caravaca et al. 2009). A recent study carried out homeostasis (Matschinsky & Wilson 2019). Another DE by Pizarro Inostroza et al. (2019) supported the strong gene of interest was OGDH, which encodes oxoglutarate association of the CSN1S1 and CSN3 genotypes with milk © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Selective sweeps in Murciano-Granadina goats 689 Figure 3 (a) Sequencing electropherograms showing the region containing the c.801C>G, (p.Cys267Trp) SNP in individuals with CC, GC and GG genotypes. (b) Cluster plots of TaqMan genotyping results obtained with the Genotyping Analysis Module implemented in the ThermoFisher Cloud computing application (Applied Biosystems). The horizontal and vertical axes correspond to alleles C and G respectively. The dots with red, green, and blue colors represent CC, CG and GG genotypes respectively. The negative control is indicated by an orange dot. Five samples with known genotypes (GG, N = 1; CC, N = 2; and GC, N = 2) assessed through Sanger sequencing were used as positive controls. (c) Manhattan plot depicting associations between coat color (41 brown and 49 black Murciano-Granadina goats) and the genotypes of marker rs669694251 (red dot) plus 1134 additional Goat SNP50 BeadChip markers mapping to goat chromosome 18. The dashed line represents the negative log10 P-value defining the threshold of significance (q ≤ 0.05) after correcting for multiple testing with a false discovery rate approach (Benjamini & Hochberg 1995). Significant SNPs are indicated with blue dots. composition traits in MG goats. In the light of these results, Despite its modest statistical significance our finding that artificial selection has targeted casein genes (q = 1.4 × 10−3), the chromosome 21 (19.88–21.86 Mb) in MG goats makes sense. Indeed, casein genotypes have selective sweep is also worth to mention because it been used since 2012 in the breeding program of the MG contains two genes that are DE in milking and dried breed as selection criteria (Delgado et al. 2017). goats i.e. PLIN1, which is involved in the synthesis of lipid The most significant SNP mapping to the chromosome 16 droplets (Itabe et al. 2017), and ANPEP, a proteolytic selective sweep (40.09–51.83 Mb, q = 1.05 × 10−4) lies enzyme (Turner 2004). Another significant chromosome close to the solute carrier family 45 member 1 (SLC45A1) 21 selective sweep (5.85–8.98 Mb, q = 0.01) encompasses gene, which encodes a glucose transporter (Vitavska & the IGF1R gene, which is also DE in lactating and dried Wieczorek 2013). In this region, we also found the β-1,3- goats (Guan et al. 2020). This locus encodes a receptor galactosyltransferase 6 (B3GALT6) gene, which mediates binding insulin-like growth factor 1, a key mediator of the transfer of galactose. Both glucose and galactose are mammary terminal end bud and ductal formation during necessary for the synthesis of lactose (Hennet 2002). development and also an essential player in mammary Besides, this region contains several metabolic loci that are gland function and maintenance (Christopoulos et al. DE in lactating and non-lactating goats (Guan et al. 2020). 2015). On chromosome 26 (26.69–31.05 Mb), we also One of these genes encodes phosphogluconate dehydroge- detected a highly significant selective sweep nase, the second dehydrogenase in the pentose phosphate (q = 4.11 × 10−8). Genes of interest mapping to this shunt (Berdis & Cook 1993). This biochemical route is used region are SCD, which encodes an enzyme catalyzing the by cells to synthesize NADPH, a molecule essential for the synthesis of monounsaturated fatty acids (Miyazaki & reductive biosynthesis of lipids (Giacomini et al. 2020). We Ntambi 2003), and NFKB2, which is a key regulator of also detected the ENO1 gene, which yields a glycolytic milk production (Connelly et al. 2010). Noteworthy, both enzyme (Kang et al. 2008), and PIK3CD, a locus involved in loci are DE in lactating and non-lactating goats (Guan the regulation of insulin signaling (Du et al. 2018). et al. 2020). © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 690 Guan et al. About the genomic architecture of pigmentation in MG In conclusion, by integrating a selection scan and other goats sources of information we have been able to identify several genes, with important roles in metabolism and lactation, Several genes with well-defined roles in pigmentation which might have been targeted by selection for dairy traits. (http://www.espcr.org/micemut/, accessed 27 February One of the most compelling cases would be that of the CSN3 2020) co-localized with the selective sweeps identified in and CSN1S2 genes, which are involved in casein synthesis. our study. The GLI3 gene, for instance, maps to the highly Moreover, we have demonstrated that several pigmentation significant chromosome 4 selective sweep (41.8–49.95 Mb) genes co-localized with selective sweeps but, according to and its inactivation leads to hypopigmentation due to a loss our data and previous results (Fontanesi et al. 2009), the of melanoblast specification (Matera et al. 2008). Another coat coloration of MG goats fundamentally depends on gene of interest is MITF, which co-localizes with the MC1R genotype. chromosome 22 (31.75 Mb) selective sweep. The MITF transcription factor regulates the expression of numerous pigmentation genes to promote melanocyte differentiation Acknowledgements (Kawakami & Fisher 2017) and its polymorphism has been This research was funded by the European Regional associated with white spotting phenotypes in horses Development Fund (FEDER)/Ministerio de Ciencia e Inno- (Hauswirth et al. 2012), and dogs (Baranowska Körberg vación – Agencia Estatal de Investigación/Project Refer- et al. 2014) as well as with piebaldism in cattle (Fontanesi ence grant: PID2019-105805RB-I00 and by the CERCA et al. 2012). A third pigmentation gene of obvious interest is Programme/Generalitat de Catalunya. We also acknowl- MC1R, which maps to a highly significant edge the support of the Spanish Ministry of Economy and (q = 7.19 × 10−6) selective sweep on chromosome 18 Competitivity for the Center of Excellence Severo Ochoa (16.64–19.46 Mb). Importantly, the melanocyte receptor 2020–2023 (CEX2019-000902-S) grant awarded to the encoded by this gene is a major regulator of the synthesis of Centre for Research in Agricultural Genomics (CRAG, melanin (Wolf Horrell et al. 2016). Bellaterra, Spain). Dailu Guan was funded by a PhD Given that our study identified several pigmentation fellowship from the China Scholarship Council (CSC). genes which might have been targeted by selection, we Marı́a Gracia Luigi-Sierra was funded with a fellowship aimed to elucidate their potential involvement in the coat Formación de Personal Investigador (BES-C-2017- coloration of MG goats. To achieve this goal, we carried out 079709) awarded by the Spanish Ministry of Economy a GWAS comprising 529 MG goats with color phenotypes. and Competitivity. The authors are indebted to the We identified a single peak surpassing the threshold of Asociación Nacional de Criadores de Caprino de Raza significance (q ≤ 0.05) on chromosome 18. This region Murciano-Granadina (CAPRIGRAN), and more specifically contains the MC1R gene, thus confirming one of the most to Miguel Garcı́a, Teresa Novoa and Manolo Delgado, for relevant results of the selection scan. Importantly, a providing the blood samples and color phenotypes of previous study based on a candidate gene approach goats analysed in the current study. We also want to reported a strong association between MC1R genotype acknowledge Drs. Juan Manuel Serradilla (Universidad de (c.801C>G, p.Cys267Trp) and coat color in 28 MG goats Córdoba), Baltasar Urrutia (Instituto Murciano de Inves- (Fontanesi et al. 2009). tigación y Desarrollo Agrario y Alimentario), and Juan Subsequent Sanger sequencing and TaqMan genotyping Carrizosa (Instituto Murciano de Investigación y Desar- experiments confirmed that the coat color of MG goats is rollo Agrario y Alimentario) for providing pictures of fully explained by MC1R c.801C>G genotype (all brown Murciano-Granadina goats. goats are CC and all black goats are either GG or GC), which is consistent with the results reported by Fontanesi et al. (2009). The G-allele (Trp) might be a dominant gain-of- Data availability statement function mutation located in the extracellular loop 3 of the The dataset used to perform GWAS and selection scans is MC1R. We hypothesize that its mechanism of action might accessible at Figshare (https://doi.org/10.6084/m9.figsha involve the disruption of a disulfide bond between Cys267 re.11999823), and MC1R sequences have been deposited and Cys275, thus altering the three-dimensional structure in the GenBank database (accession codes: MT186757– and activity of the MC1R protein towards an increased MT186778). production of eumelanin (Wolf Horrell et al. 2016), but measuring the activity of the Cys267 and Trp267 MC1R alleles in cultured melanocytes would be essential to References confirm this. Dominant mutations increasing eumelanin Bâlteanu V.A., Cardoso T.F., Amills M., Luigi-Sierra M.G., Egerszegi synthesis have been previously reported in pigs, cattle, and I., Anton I. & Zsolnai A. (2021) Red and blond Mangalitza pigs sheep (Klungland et al. 1995; Kijas et al. 1998; Våge et al. display a signature of divergent directional selection in the 1999). SLC45A2 gene. Animal Genetics 52, 66–77. © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Selective sweeps in Murciano-Granadina goats 691 Baranowska Körberg I., Sundström E., Meadows J.R.S. et al. (2014) Fontanesi L., Beretti F., Riggio V., Dall’Olio S., González E.G., A simple repeat polymorphism in the MITF-M promoter is a key Finocchiaro R., Davoli R., Russo V. & Portolano B. (2009) regulator of white spotting in dogs. PLoS One 9, e104363. Missense and nonsense mutations in melanocortin 1 receptor Benjamini Y. & Hochberg Y. (1995) Controlling the false discovery (MC1R) gene of different goat breeds: association with red and rate – a practical and powerful approach to multiple testing. black coat colour phenotypes but with unexpected evidences. Journal of the Royal Statistical Society Series B-Statistical Method- BMC Genetics 10, 47. ology 57, 289–300. Fontanesi L., Scotti E. & Russo V. (2012) Haplotype variability in Berdis A.J. & Cook P.F. (1993) Overall kinetic mechanism of 6- the bovine MITF gene and association with piebaldism in phosphogluconate dehydrogenase from Candida utilis. Biochem- Holstein and Simmental cattle breeds. Animal Genetics 43, 250–6. istry 32, 2036–40. Gautier M., Klassmann A. & Vitalis R. (2017) rehh 2.0: a Bertolini F., Servin B., Talenti A. et al. (2018) Signatures of reimplementation of the R package rehh to detect positive selection and environmental adaptation across the goat genome selection from haplotype structure. Molecular Ecology Resources post-domestication. Genetics Selection Evolution 50, 57. 17, 78–90. Bickhart D.M., Rosen B.D., Koren S. et al. (2017) Single-molecule Gautier M. & Vitalis R. (2012) rehh: an R package to detect sequencing and chromatin conformation capture enable de novo footprints of selection in genome-wide SNP data from haplotype reference assembly of the domestic goat genome. Nature Genetics structure. Bioinformatics 28, 1176–7. 49, 643–50. Giacomini I., Ragazzi E., Pasut G. & Montopoli M. (2020) The Brito L.F., Kijas J.W., Ventura R.V., Sargolzaei M., Porto-Neto L.R., pentose phosphate pathway and its involvement in cisplatin Cánovas A., Feng Z., Jafarikia M. & Schenkel F.S. (2017) Genetic resistance. International Journal of Molecular Sciences 21, 937. diversity and signatures of selection in various goat breeds Guan D., Landi V., Luigi-Sierra M.G. et al. (2020) Analyzing the revealed by genome-wide SNP markers. BMC Genomics 18, 229. genomic and transcriptomic architecture of milk traits in Caravaca F., Amills M., Jordana J. et al. (2008) Effect of αs1-casein Murciano-Granadina goats. Journal of Animal Science and Biotech- (CSN1S1) genotype on milk CSN1S1 content in Malagueña and nology 11, 35. Murciano-Granadina goats. Journal of Dairy Research 75, 481–4. E G.-X., Zhao Y.-J. & Huang Y.-F. (2019) Selection signatures of Caravaca F., Ares J.L., Carrizosa J. et al. (2011) Effects of αs1-casein litter size in Dazu black goats based on a whole genome (CSN1S1) and κ-casein (CSN3) genotypes on milk coagulation sequencing mixed pools strategy. Molecular Biology Reports 46, properties in Murciano-Granadina goats. Journal of Dairy Research 5517–23. 78, 32–7. Harris A.M., Garud N.R. & DeGiorgio M. (2018) Detection and Caravaca F., Carrizosa J., Urrutia B., Baena F., Jordana J., Amills M., classification of hard and soft sweeps from unphased genotypes Badaoui B., Sánchez A., Angiolillo A. & Serradilla J.M. (2009) by multilocus genotype identity. Genetics 210, 1429–52. Short communication: effect of αS1-casein (CSN1S1) and κ-casein Hauswirth R., Haase B., Blatter M. et al. (2012) Mutations in MITF (CSN3) genotypes on milk composition in Murciano-Granadina and PAX3 cause “Splashed White” and other white spotting goats. Journal of Dairy Science 92, 2960–4. phenotypes in horses. PLoS Genetics 8, e1002653. Christopoulos P.F., Msaouel P. & Koutsilieris M. (2015) The role of Hennet T. (2002) The galactosyltransferase family. Cellular and the insulin-like growth factor-1 system in breast cancer. Molec- Molecular Life Sciences CMLS 59, 1081–95. ular Cancer 14, 43. Huang D.W., Sherman B.T. & Lempicki R.A. (2008) Systematic and Connelly L., Barham W., Pigg R., Saint-Jean L., Sherrill T., Cheng integrative analysis of large gene lists using DAVID bioinformat- D.-S., Chodosh L.A., Blackwell T.S. & Yull F.E. (2010) Activation ics resources. Nature Protocols 4, 44–57. of nuclear factor kappa B in mammary epithelium promotes milk Huang D.W., Sherman B.T. & Lempicki R.A. (2009) Bioinformatics loss during mammary development and infection. Journal of enrichment tools: paths toward the comprehensive functional Cellular Physiology 222, 73–81. analysis of large gene lists. Nucleic Acids Research 37, 1–13. de Oliveira P.SN., Cesar A.SM., do Nascimento M.L. et al. (2014) Itabe H., Yamaguchi T., Nimura S. & Sasabe N. (2017) Perilipins: a Identification of genomic regions associated with feed efficiency diversity of intracellular lipid droplet proteins. Lipids in Health and in Nelore cattle. BMC Genetics 15, 100. Disease 16, 83. Delgado J.V., Landi V., Barba C.J., Fernández J., Gómez M.M., Kang H.J., Jung S.K., Kim S.J. & Chung S.J. (2008) Structure of Camacho M.E., Martı́nez M.A., Navas F.J. & León J.M. (2017) human α-enolase (hENO1), a multifunctional glycolytic enzyme. Murciano-Granadina Goat: A Spanish local breed ready for the Acta Crystallographica Section D: Biological Crystallography 64, challenges of the twenty-first century. In: Sustainable Goat 651–7. Production in Adverse Environments: Volume II: Local Goat Breeds Kawakami A. & Fisher D.E. (2017) The master role of (Ed. by J. Simões & C. Gutiérrez), pp. 205–19. Springer microphthalmia-associated transcription factor in melanocyte International Publishing, Cham, Switzerland. and melanoma biology. Laboratory Investigation 97, 649–56. Du X., Li X., Chen L. et al. (2018) Hepatic miR-125b inhibits insulin Kijas JMH, Wales R., Törnsten A., Chardon P., Moller M. & signaling pathway by targeting PIK3CD. Journal of Cellular Andersson L. (1998) Melanocortin receptor 1 (MC1R) mutations Physiology 233, 6052–66. and coat color in pigs. Genetics 150, 1177–85. Flint D.J., Tonner E., Knight C.H., Whitelaw CBA, Webster J., Barber Kim E.S., Elbeltagy A.R., Aboul-Naga A.M., Rischkowsky B., Sayre M. & Allan G. (2001) Control of mammary involution by insulin- B., Mwacharo J.M. & Rothschild M.F. (2016) Multiple genomic like growth factor binding proteins: role of prolactin. Livestock signatures of selection in goats and sheep indigenous to a hot arid Production Science 70, 115–20. environment. Heredity 116, 255–64. © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 692 Guan et al. Kim H.-S. (2013) Role of insulin-like growth factor binding protein- content and quality in Murciano-Granadina goats. Journal of 3 in glucose and lipid metabolism. Annals of Pediatric Endocrinol- Animal Breeding and Genetics 137, 407–22. ogy & Metabolism 18, 9–12. Purcell S., Neale B., Todd-Brown K. et al. (2007) PLINK: a tool set Klungland H., Vage D.I., Gomez-Raya L., Adalsteinsson S. & Lien S. for whole-genome association and population-based linkage (1995) The role of melanocyte-stimulating hormone (MSH) analyses. The American Journal of Human Genetics 81, 559–75. receptor in bovine coat color determination. Mammalian Genome Quinlan A.R. & Hall I.M. (2010) BEDTools: a flexible suite of utilities 6, 636–9. for comparing genomic features. Bioinformatics 26, 841–2. Kumar S., Stecher G., Li M., Knyaz C. & Tamura K. (2018) MEGA X: Saravanan K.A., Panigrahi M., Kumar H., Bhushan B., Dutt T. & Molecular Evolutionary Genetics Analysis across computing Mishra B.P. (2020) Selection signatures in livestock genome: a platforms. Molecular Biology and Evolution 35, 1547–9. review of concepts, approaches and applications. Livestock Science Lai F.-N., Zhai H.-L., Cheng M. et al. (2016) Whole-genome 241, 104257. scanning for the litter size trait associated genes and SNPs under Song S., Yao N., Yang M. et al. (2016) Exome sequencing reveals selection in dairy goat (Capra hircus). Scientific Reports 6, 38096. genetic differentiation due to high-altitude adaptation in the Li X., Su R., Wan W. et al. (2017) Identification of selection signals Tibetan cashmere goat (Capra hircus). BMC Genomics 17, 122. by large-scale whole-genome resequencing of cashmere goats. Tosser-Klopp G., Bardou P., Bouchez O. et al. (2014) Design and Scientific Reports 7, 15142. characterization of a 52K SNP Chip for goats. PLoS One 9, Liao Y., Smyth G.K. & Shi W. (2014) featureCounts: an efficient e86227. general purpose program for assigning sequence reads to Turner A.J. (2004) Membrane alanyl aminopeptidase. In: Handbook genomic features. Bioinformatics 30, 923–30. of Proteolytic Enzymes (2nd edn) (Ed by A.J. Barrett, N.D. Love M.I., Huber W. & Anders S. (2014) Moderated estimation of Rawlings & J.F. Woessner), pp. 289–94. Academic Press, London, fold change and dispersion for RNA-seq data with DESeq2. UK. Genome Biology 15, 550. Våge D.I., Klungland H., Lu D. & Cone R.D. (1999) Molecular and Martin P., Palhiere I., Maroteau C. et al. (2017) A genome scan for pharmacological characterization of dominant black coat color in milk production traits in dairy goats reveals two new mutations sheep. Mammalian Genome 10, 39–43. in Dgat1 reducing milk fat content. Scientific Reports 7, 1872. Vitavska O. & Wieczorek H. (2013) The SLC45 gene family of Martı́nez-Reyes I. & Chandel N.S. (2020) Mitochondrial TCA cycle putative sugar transporters. Molecular Aspects of Medicine 34, metabolites control physiology and disease. Nature Communica- 655–60. tions 11, 102. Vitti J.J., Grossman S.R. & Sabeti P.C. (2013) Detecting natural Matera I., Watkins-Chow D.E., Loftus S.K., Hou L., Incao A., Silver selection in genomic data. Annual Review of Genetics 47, 97–120. D.L., Rivas C., Elliott E.C., Baxter L.L. & Pavan W.J. (2008) A Voight B.F., Kudaravalli S., Wen X. & Pritchard J.K. (2006) A map sensitized mutagenesis screen identifies Gli3 as a modifier of of recent positive selection in the human genome. PLoS Biology 4, Sox10 neurocristopathy. Human Molecular Genetics 17, 2118–31. e72. Matschinsky F.M. & Wilson D.F. (2019) The central role of Wolf Horrell E.M., Boulanger M.C. & D’Orazio J.A. (2016) glucokinase in glucose homeostasis: a perspective 50 years after Melanocortin 1 receptor: structure, function, and regulation. demonstrating the presence of the enzyme in islets of Langer- Frontiers in Genetics 7, 95. hans. Frontiers in Physiology 10, 148. Xu S.-S., Ren X., Yang G.-L. et al. (2017) Genome-wide association Miyazaki M. & Ntambi J.M. (2003) Role of stearoyl-coenzyme A analysis identifies the genetic basis of fat deposition in the tails of desaturase in lipid metabolism. Prostaglandins, Leukotrienes and sheep (Ovis aries). Animal Genetics 48, 560–9. Essential Fatty Acids 68, 113–21. Zhang B., Chang L., Lan X. et al. (2018) Genome-wide definition of Mychaleckyj J.C., Nayak U., Colgate E.R. et al. (2018) Multiplex selective sweeps reveals molecular evidence of trait-driven genomewide association analysis of breast milk fatty acid domestication among elite goat (Capra species) breeds for the composition extends the phenotypic association and potential production of dairy, cashmere, and meat. GigaScience 7, giy105. selection of FADS1 variants to arachidonic acid, a critical infant Zhou X. & Stephens M. (2012) Genome-wide efficient mixed-model micronutrient. Journal of Medical Genetics 55, 459–68. analysis for association studies. Nature Genetics 44, 821–4. Neuenschwander S., Schwartz A., Wood T.L., Roberts C.T. Jr, Hennighausen L. & LeRoith D. (1996) Involution of the lactating Supporting information mammary gland is inhibited by the IGF system in a transgenic mouse model. The Journal of Clinical Investigation 97, 2225–32. Additional supporting information may be found online in Oget C., Servin B. & Palhière I. (2018) Genetic diversity analysis of the Supporting Information section at the end of the article. French goat populations reveals selective sweeps involved in their Figure S1 Description of the experimental design. differentiation. Animal Genetics 50, 54–63. Table S1 Primers used in the amplification of the coding Onzima R.B., Upadhyay M.R., Mukiibi R., Kanis E., Groenen M.A.M. region of the melanocortin 1 receptor gene. & Crooijmans R.P.M.A. (2018) Genome-wide population struc- Table S2 Custom TaqMan probes used to genotype the ture and admixture analysis reveals weak differentiation among Ugandan goat breeds. Animal Genetics 49, 59–70. rs669694251 polymorphism. Pizarro Inostroza M.G., Landi V., Navas González F.J., León Jurado Table S3 List of selective sweeps detected with the iHS J.M., Martı́nez Martı́nez M.D.A., Fernández Á.J. & Delgado statistic in Murciano-Granadina goats. Bermejo J.V. (2019) Non-parametric association analysis of Table S4 List of enriched pathways/terms using the set of additive and dominance effects of casein complex SNPs on milk genes co-localizing with putative selective sweeps. © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693 13652052, 2021, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/age.13113 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [21/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Selective sweeps in Murciano-Granadina goats 693 Table S5 List of genes mapping to selective sweeps and Table S6 Genotypes of the rs657434682 and rs669694251 showing differential expression in lactating (T1 and T2) and SNP markers assessed through the visual inspection of non-lactating (T3) goats. sequencing electropherograms. © 2021 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 52, 683–693