This is the accepted version of the article: Luigi Sierra, Maria Gracia; Landi, Vicenzo; Guan, Dailu; [et al.]. «A genome- wide association analysis for body, udder, and leg conformation traits recorded in Murciano-Granadina goats». Journal of Dairy Science, Vol. 103, Issue 12 (December 2020), p. 11605-11617. DOI 10.3168/jds.2020-18461 This version is available at https://ddd.uab.cat/record/238716 under the terms of the license Interpretive summary: The genetic basis of morphological traits has not been investigated in depth in goats despite its great impact on productive performance. Herewith, we have genotyped 825 Murciano-Granadina goats with records for 17 morphological traits with a chip containing 54,241 markers. Statistical analyses have made possible to identify two genomic regions displaying highly significant associations with medial suspensory ligament (chromosomes 16 and 28). This is a first step towards identifying the genetic factors that influence morphological traits with the goal of incorporating such information to breeding plans. RUNNING HEAD: GWAS FOR MORPHOLOGY TRAITS IN GOATS A genome-wide association analysis for body, udder and leg conformation traits recorded in Murciano-Granadina goats Maria Gracia Luigi-Sierra1, Vincenzo Landi2,3, Dailu Guan1, Juan Vicente Delgado2, Anna Castelló1,4, Betlem Cabrera1,4, Emilio Mármol-Sánchez1, Javier Fernández Alvarez5, Mayra Gómez-Carpio2, Amparo Martínez2, Xavier Such6, Jordi Jordana4, Marcel Amills1,4, * 1 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Campus Universitat Autònoma de Barcelona, Bellaterra, 08193, Spain. 2 Departamento de Genética, Universidad de Córdoba, Córdoba 14071, Spain. 3 Department of Veterinary Medicine, University of Bari "Aldo Moro", SP. 62 per Casamassima km. 3, 70010 Valenzano (BA), Italy. 4 Departament de Ciència Animal i dels Aliments, Facultat de Veterinària, Universitat Autònoma de Barcelona, Bellaterra 08193, Spain, 5 Asociación Nacional de Criadores de Caprino de Raza Murciano-Granadina (CAPRIGRAN), 18340 Granada, Spain. 6 Group of Research in Ruminants (G2R), Department of Animal and Food Science, Universitat Autònoma de Barcelona (UAB), Bellaterra, Barcelona, Spain 1 Corresponding author: marcel.amills@uab.cat 2 ABSTRACT Morphological traits are of great importance to dairy goat production given their incidence on phenotypes of economic interest. However, their genomic architecture has not yet been extensively characterized. Herewith we aimed to identify genomic regions associated with 17 body, udder and leg conformation traits recorded in 825 Murciano- Granadina goats. We have gernotyped this resource population with the Illumina GoatSNP50 Beadchip (Illumina inc, San Diego, CA) and we have performed genome- wide association analyses with the GEMMA software. By doing so, we have found two genome-wide significant associations between markers rs268273468 (CHI 16:69617700) and rs268249346 (CHI 28:18321523) and medial suspensory ligament. In contrast, we have not detected any genome-wide significant associations for body and leg traits. Moreover, we have found 12, 19 and 7 chromosome-wide significant associations for udder, body, and leg traits, respectively. Comparison of our data with previous studies revealed a low level of positional concordance between regions associated with morphological traits. Besides technical factors, this observation could be due to the existence of a substantial level of genetic heterogeneity amongst breeds and/or to the strong polygenic background of morphological traits, a circumstance that makes it difficult to detect genetic factors with small phenotypic effects. Key Words. GWAS, goat, Murciano-Granadina, morphological traits. 3 INTRODUCTION Since its establishment in 1975 (Delgado et al., 2018), Murciano-Granadina has become the most important dairy goat breed in Spain, reaching a census of 112,417 heads in 2019 (https://www.mapa.gob.es/es/ganaderia/temas/zootecnia/razas-ganaderas/razas). Murciano-granadina shows a great capacity to adapt to harsh environments, yielding 530 kg of milk per lactation (250 days) with fat and protein contents of 5.6% and 3.6%, respectively (https://www.mapa.gob.es/es/ganaderia/temas/zootecnia/razas- ganaderas/razas). The genetic improvement program of the Murciano-Granadina breed is mostly focused on augmenting milk yield and quality as well as on optimizing body and dairy conformation traits. The association of Murciano-Granadina breeders (Caprigran) records systematically information about milk yield and composition and they also perform the linear scoring of 17 morphological traits (Delgado et al., 2018). The inclusion of morphological traits as selection criteria is motivated by their relationship with mammary health and longevity (Shelton, 1978; Manfredi et al., 2001; Montaldo and Manfredi, 2002). Studies performed in cattle (Seykora and McDaniel, 1985; Rogers et al., 1991; Boettcher et al., 1998; Rupp and Boichard, 1999; Miglior et al., 2017) support a relationship between udder morphology and health. Indeed, higher, non-pendulous and more tightly attached udders are less susceptible to mastitis (Seykora and McDaniel, 1985; Rupp and Boichard, 2003). Moreover, flat, disk, or inverted teat ends are associated with an increased risk of suffering mastitis, while funnel shaped teats seem to be less prone to this disease (Seykora and McDaniel, 1985). There are also numerous reports that associate body traits with fertility and longevity (Bastin and Gengler 2013; Miglior et al., 2017). For instance, feet and leg conformation traits could be considered as indicators of claw health which, after reproduction and 4 mastitis, is one of the main determinants for culling animals (Egger-Danner et al., 2014). Far less studies about the genetics of morphological traits have been performed in goats. According to Rupp et al. (2011), udder floor position showed negative genetic correlations with somatic cell score both in Alpine (rg = -0.24) and Saanen goats (rg = - 0.19). The same authors demonstrated that in the Saanen breed, somatic cell score was genetically correlated with teat length (rg = 0.29), teat width (rg = 0.34), and teat form (rg = -0.27). These results suggest that a reduction in somatic cell count can be achieved by selection, while still improving milk production and udder type and teat traits. In Tinerfeña goats, moderate to high and mostly positive phenotypic correlations have been detected between udder traits and milk yield (Capote et al., 2006). Manfredi et al. (2001) found a high heritability (h2 > 0.4) for thorax perimeter in Alpine and Saanen goats, while other body traits such as rump angle, feet angle and hock distance showed low heritabilities (h2 = 0.03-0.16). In contrast, heritabilities for udder and teat scores reached values around 0.3 for most traits, with teat angle displaying the lowest value (h2 = 0.15 in Saanen). Genetic correlations among teat dimension traits and between udder floor and rear udder attachment (rg > 0.7 in Alpine and Saanen) were generally high, while the majority of genetic correlations between body and udder scores reached values below 0.3 (Manfredi et al., 2001). In a more recent study, McLaren et al. (2016) described low to moderate heritabilities (from 0.02 to 0.38) for conformation traits recorded in mixed-breed dairy goats: while udder and teat traits displayed the highest heritabilities (h2 ~ 0.28, from 0.15 to 0.38), feet and leg traits showed lower values (h2 ~ 0.13, from 0.02 to 0.25). Although the majority of the correlations estimated between milk yield and udder and teat traits were 5 negative, their magnitude and sign fluctuated across the first lactation. For instance, the estimates of the genetic correlation between udder furrow and milk yield went from −0.42 to 0.18 depending on the time point of the 1st lactation at which they were calculated (McLaren et al., 2016). Castañeda-Bustos et al. (2017) described a high genetic correlation between the productive life of dairy goats (i.e. total days in production until 72 months of age) and final score (appraisal of the general conformation of the animal), fore udder attachment and rump width. These findings demonstrate that taking into account conformation and udder traits in selection schemes could be expected to increase productivity without compromising the fitness of the animals (Castañeda-Bustos et al., 2017). Few investigations have been carried out to identify genomic regions associated with body conformation and udder traits in goats. Mucha et al. (2018) performed a genome- wide association study (GWAS) for morphology traits in mixed-breed dairy goats and found a significant association between a region on chromosome 19 and udder attachment, udder depth and front legs morphology. Moreover, Martin et al. (2018) detected 37 genome-wide significant QTL for type and somatic cell count phenotypes with linkage analyses, while a much larger number of QTL were identified by association mapping. These authors concluded that the inheritance of body and udder conformation traits is markedly polygenic and that genetic determinants are often breed- specific. In the current work, we aimed to identify genomic regions associated with the phenotypic variation of 17 morphological traits in Murciano-Granadina goats. MATERIALS AND METHODS 6 Phenotypic Recording A total of 825 female goats distributed in 13 farms, with an average herd size of 500 individuals, were scored for 17 morphological traits included in the breeding program of the Murciano-Granadina goat breed. The scoring is performed only once in the lifetime of the animal. The majority of the goats comprised in our study were scored during their first lactation, although there is a group of 89 animals that were scored between the 2nd and 6th lactation. All traits were scored by the same specialist, using a personal digital device with specific software (Kalifadroid app) for carrying out scoring tasks. The following phenotypes were evaluated, with linear scores ranging from 1 to 9 according to the measures established in Sánchez-Rodríguez (2012): (1) Udder Traits Fore Udder Attachment (FUA). Corresponds to the angle formed by the line of the udder insertion and the abdominal wall. Scores 1 and 9 correspond to angles of 45º and 120º, respectively. Rear Udder Height (RUH). This variable is scored by measuring the distance between the bottom of the vulva and the top of the secretory tissue of the mammary gland. Scores 1 and 9 correspond to distances of 11 cm and 3 cm, respectively. Udder Depth (UD). It is the distance from the lowest part of the udder floor to the hock joint (tibiotarsal joint). A linear score of 1 corresponds to an udder with its deepest part 10 cm over the hock, while a score of 9 would define an udder with its deepest part 10 cm down the hock joint 7 Medial Suspensory Ligament (MSL). Is the depth of the udder cleft measured at the base of the rear udder. Scores 1 and 9 correspond to 1 cm and 9 cm (or more) deep udder clefts, respectively. Udder Width (UW). Is measured at the crease where the udder meets the leg. Scores of 1 and 9 correspond to measurements of 3 cm and 11 cm (or more), respectively. Teat Diameter (TD). It is the diameter of the teat at its base, when it meets the udder. Score 1 = diameter of 0.5 cm and score 9 = diameter of 4.5 cm or more. Teat Placement (TP). It defines the position of the teats on the udder half. Teats located on the outside third of the udder half are scored with 1, whilst teats located very close to the medial suspensory ligament, that almost touch each other, are scored with 9. Teats on the center of the udder half, with an intermediate placement, are considered as desirable and they are scored with 5. (2) Body conformation traits Height (HT). It measures the distance from the level ground to the top of the withers. Goats with 62 cm or less receive a score of 1, while those over 78 cm are scored with 9. Chest Width (CW). Measured from the inside surface of the chest between the top of the front legs. Score 1 = 15 cm or less. Score 9 = 23 cm or more. Body Depth (BD). Is the distance between the top of the spine and bottom of the body at the beginning of the last rib. Score 1 = low depth, if the beginning of the last rib 8 is located above the elbow joint, and score 9 = high depth, if the beginning of the last rib is located below the elbow joint. Rump Width (RW). Distance between the most posterior points of the pin bones (ischial tuberosities). Score 1 = 13 cm (or less) and score 9 = 21 cm (or more). Rump Angle (RA). Angle between the hook (coxal tuberosity) and pin (ischial tuberosity) bones. Scores 1 and 9 correspond to angles of 55º and 31º (approximately), respectively. Angularity or Dairyness (ANG). Is the angle and openness of the ribs. Score 1 defines an animal extremely coarse for this trait and score 9 is assigned to goats that are very angular. Bone Quality (BQ). It is an appraisal of the thickness and width of the bone structure, assessed by examining the rear leg from the rear and from the side. Score 1 corresponds to goats with thick and round bones, while 9 is assigned to goats with flat and sharp bones. (3) Feet structure Rear Legs Rear View (RLR). Direction of rear feet when viewed from the rear. Score 1: extreme toe-out feet and score 9: parallel feet. Rear Legs Side View (RLS). Curvature of the hock viewed laterally. Score 1: straight legs and score 9: very curved legs. Mobility (MOB). Evaluates the locomotion patterns, including the length and direction of the step. Score 1: bad locomotion, with severe abduction and short steps, and score 9: harmonic, long and uniform locomotion. 9 Isolation of Genomic DNA from Blood and Genotyping with the Goat SNP50 Beadchip Blood samples from the 825 Murciano-Granadina goats with morphology records were collected in vacuum tubes coated with EDTA K3 anticoagulant and stored at -20 ºC until processing. The purification of genomic DNA was achieved with a modified salting-out procedure (Miller et al., 1988). To do so, we combined 3 mL of whole blood plus 4 volumes of Red Cell Lysis Solution (Tris-HCl 10mM, pH = 6.5; EDTA 2 mM; Tween 20 1%) and then this mixture was centrifuged at a speed of 2000 g. The resulting cell pellet was resuspended in 3 mL of lysis buffer (Tris-HCl 200mM, pH = 8, EDTA 30 mM, SDS 1%; NACl 250 mM) and 100 µl proteinase K (20 mg/mL) and incubated during 3 hours at 55° C. The lysate was chilled, and 1 mL of ammonium acetate 10 M was added. Next, a centrifugation step at 2000 g during 10 minutes was carried out, and the supernatant (~4 mL) was transferred to a new tube with 3 mL of isopropanol 96%. Subsequently, samples were centrifuged at 2000 g for 3 minutes. The resulting DNA pellet was washed with 3 mL of ethanol 70% and an additional centrifugation step at 2000 g for 1 minute was performed. The DNA pellet was left at room temperature until it dried, and it was subsequently resuspended in 1 mL of TE buffer (Tris-HCl 10 mM, EDTA 1 mM, pH = 8). All 825 goats were genotyped with the Illumina Goat SNP50 BeadChip (Illumina Inc., San Diego, CA), which contains 54,241 SNPs, following the instructions of the manufacturer. The genomic location of the SNPs was obtained by using the goat ARS1 genome (Bickhart et al., 2017) as reference, and the position and the name of each SNP was updated using the software PLINK v 1.9 (Chang et al., 2015). The genotypic information was filtered using PLINK v 1.9 (Chang et al., 2015). Only individuals with less than 5% of missing genotypes were taken into consideration. With regard to SNPs, only those meeting the following requirements were used in the 10 GWAS: (1) Mapping to autosomes, (2) Displaying a minor allele frequency of 0.05 or higher, (3) Not deviating significantly (P < 0.001) from the Hardy-Weinberg expectation, and (4) With a genotype call rate over 90%. After applying these filtering criteria, 47,880 SNPs and 811 animals were selected to perform subsequent analyses. Population structure was assessed with a principal component analysis (PCA) approach implemented in PLINK v1.9 (Chang et al., 2015). The visualization of the PCA results was based on the first two components of the PCA. Statistical Analyses We calculated summary statistics for each one of the morphological traits using R (R core Team, 2017) as well as Pearson correlations (rp) between conformation traits and milk composition and yield records, i.e. Total milk yield (RMY), milk yield at 210 days (MY210), milk yield at 240 days (MY240), milk yield at 305 days (MY305), somatic cell count (SCC), fat percentage (FP), protein percentage (PP) and lactose percentage (LP). The correlations were estimated using R software (R Core Team, 2017). Heatmap plots were constructed with R software (R Core Team, 2017) to visualize the correlation matrix and the P-values of each correlation. The software GEMMA (Zhou and Stephens, 2012) was used to carry out the GWAS. This method corrects population structure by considering the relatedness matrix, which is built by taking into account all genome-wide SNPs as a random effect. Morphological phenotypes were rank-based transformed by using the package GenAbel from R (Aulchenko et al., 2007) because we assessed, with our data, that this transformation yields residuals that are normally distributed. A univariate linear mixed model was fit for each trait as follows, y = Wα + xβ + u + ϵ 11 where y is a vector of corrected scores for morphological traits recorded in 811 individuals; W = (w1, · · · , wc) is an n×c matrix of three fixed effects (farm, with 13 levels; year of birth, with 10 levels; number of lactations, with 6 levels) and one covariate (days producing milk); α is a c-vector of the corresponding fixed effects including the intercept; x is a n-vector of marker genotypes; β is the effect size of the marker (allele substitution effect); u is a n-vector of random individual genetic effects with a normal distribution u ∼ N (0, λ τ− 1 K), being τ− 1 the variance of the residual error, λ the ratio between the two variance components and K the relatedness matrix derived from SNP genotypes. Finally, ϵ is a n-vector of errors. A false discovery rate (FDR) approach was applied to correct for multiple testing, setting the significance level to a q-value of 0.05 (Benjamini and Hochberg, 1995). Graphical visualization of the results of the GWAS was achieved by using the R software (R core Team, 2017). The proportion of the phenotypic variance explained by the significant SNPs (PVE) was estimated using the formula reported by Shim et al. (2015), 2β̂2 MAF (1 − MAF) 𝑃𝑉𝐸 = 2β̂2 MAF (1 − MAF) + (𝑠𝑒(β̂))2 2N MAF (1 − MAF) Where β is the effect size of the SNP variant estimated from the association analysis; MAF is the minor allele frequency of the SNP; se is the standard error and N is the sample size. Lambda inflation factors were calculated with the median method (1 d.f.) implemented in GenABEL (Aulchenko et al. 2007), while Q-Q plots were built with the gg_qqplot() function (https://www.rdocumentation.org/packages/lindia/versions/0.9 /topics/gg_qqplot). 12 In order to retrieve candidate genes mapping close to significant SNPs we considered an interval of ± 1 Mb based on data previously reported by Guan et al., (2020) for the same population. Genes mapping within these defined boundaries were listed by using the Biomart tool from Ensemble (Kinsella et al., 2011), and subsequently analyzed with Uniprot (Bateman, 2019) and David Bioinformatic Resources (Huang et al., 2009) to annotate their function. RESULTS AND DISCUSSION Population Structure and Analysis of Morphological Traits The first two components C1 and C2 of the principal component analysis accounted for 42.64% and 21.76 % of the genetic variance, respectively (Supplementary Figure 1). A number of samples grouped according to their farm of origin, but in general we did not find an obvious within-population substructure in this sample of Murciano- Granadina goats. Descriptive statistics of the raw conformation scores are reported in Supplementary Table 1 and Supplementary Figures 2 and 3. The estimated phenotypic correlations and their significances are depicted 1n Supplementary Table 2 and Figure 2. We have classified phenotypic correlations as low (rP < 0.2), moderate (rP = 0.2-0.4) or high (rP > 0.4). Phenotypic correlations between udder traits were in general low and positive, with the exception of the correlation between medial suspensory ligament and udder depth that was moderate (rP = 0.369, P-value < 0.001). In the study of McLaren et al. (2016) genetic correlations between udder traits ranged from 0.12 to 0.77, while those between teat traits were in the range of -0.10 to 0.69. 13 Udder and teat traits were positively correlated, a result coincident with our findings, and also with those reported by Manfredi et al. (2001). We found low and positive correlations between somatic cell score and medial suspensory ligament (rp = 0.14, P- value < 0.001), UD (rp = 0.189, P-value < 0.001) and TD (rp = 0.131, P-value < 0.001). Udder depth (rG = 0.10) and teat size (rG = 0.29) also showed positive genetic correlations with somatic cell score in Latxa sheep (Legarra and Ugarte, 2005). It is well known that pendulous udders are closer to the floor, exposing the mammary gland to direct contact with fecal and other environmental contaminants (Pugh and Baird, 2012). Moreover, teat diameter is highly correlated with the diameter of the internal cistern (Guarín et al., 2017), and it is expected that sphincters in animals with wider cisterns do not close completely, leaving a channel open for pathogens, thus increasing the risk of suffering mastitis (Seykora and McDaniel, 1985). In contrast, udder width (rP = -0.109, P-value < 0.001) and teat placement (rP = -0.126, P-value < 0.001) were negatively correlated with somatic cell score, but these two values were low and in other studies these two parameters do not show strong genetic correlations with somatic cell score (Legarra and Ugarte, 2005; Pérez-Cabal et al., 2013). Notably, udder width (rP = 0.20- 0.21, P-values < 0.001) and udder depth (rP = 0.29-0.32, P-values < 0.001) showed moderate positive correlations with milk yield. Pérez-Cabal et al. (2013) found moderate to high and positive phenotypic correlations between milk yield and udder width (rP = 0.29) and udder depth (rP = 0.47) in Spanish Assaf sheep, and Legarra and Ugarte (2005) reported a strong positive genetic correlation between milk yield and udder depth (rG = 0.43). In the end, it is reasonable to infer that goats with wider and deeper udders produce more milk. Phenotypic correlations between body conformation traits were generally high and positive. For instance, height was correlated with chest width (rP = 0.646, P-value < 14 0.001), rump width (rP = 0.673, P-value < 0.001) and angularity (rP = 0.424, P-value < 0.001), and chest width also showed positive correlations with body depth (rP = 0.415, P-value < 0.001), rump width (rP = 0.744, P-value < 0.001) and angularity (rP = 0.682, P-value < 0.001). Zujovic et al. (2011) also observed high phenotypic correlations between body traits measured in Balkan goats, and a similar trend was observed by Chacón et al. (2011) in Cuban goats. In general, taller goats are also bigger and have a wider chest and rump. We also observed moderate positive correlations between chest width and body depth and milk yield (rP = 0.2-0.3, P-values < 0.001) and between rump width and milk yield (rP =0.15-0.2, P-values < 0.001). In cows it has been reported that animals with a wider rump and increased space across their hooks and pins can accommodate a higher and wider udder Campbell and Marshall (2016). These results suggest that the increment of the udder capacity leads to a higher milk yield while at the same time it decreases the percentages of solid milk components (% lactose, % fat and % protein). There are also indications that cows with a short and round body often lack dairy character and udder capacity (Campbell and Marshall, 2016), supporting the observed positive correlation between angularity and milk yield in our Murciano- Granadina population. With regard to leg traits, mobility displayed a positive phenotypic correlation with rear legs rear view (rP = 0.425, P-value < 0.001) and RLS (rP = 0.231, P-value < 0.001). Interestingly rear legs rear view and rear legs side view showed low positive and low negative correlations with milk yield, respectively. McLaren et al. (2016) observed a correlation of 0.33 between back legs and milk yield at 305 days, while de la Fuente et al. (2011) estimated a small genetic correlation of −0.09 between the back legs and milk yield of Churra ewes. Our interpretation is that leg morphology could be associated with 15 the predisposition of goats to suffer lameness, a pathology that results in a decreased milk production and often in the culling of the affected animal (Archer et al., 2010). Identification of Genetic Determinants for Udder Traits After performing the GWAS analysis we detected two genome-wide significant associations for the trait medial suspensory ligament (Table 1). We also found chromosome-wide associations for the traits medial suspensory ligament (2 SNPs), udder width (1 SNPs), udder depth (3 SNPs), teat placement (4 SNPs) and teat diameter (2 SNPs), as shown in Table 1 and Supplementary Figures 4 (Manhattan plots) and 7 (Q-Q plots). No significant SNPs were found for fore udder attachment and rear udder height. The rs268273468 (CHI 16: 69617700) marker, which was significantly associated to medial suspensory ligament at the genome-wide level of significance (Table 1 and Figure 1), is located less than 1 Mb away from the lysophosphatidylglycerol acyltransferase 1 (LPGAT1) gene (Supplementary Table 3). This gene encodes an enzyme involved in the conversion of lysophosphatidylglycerol into phosphatidylglycerol, a membrane phospholipid that is a key precursor in the biosynthesis of cardiolipin (Yang et al., 2004). Interestingly, cardiolipin is located in the inner mitochondrial membrane and plays a fundamental role in maintaining mitochondrial membrane stability and dynamics, as well as regulating apoptosis (Paradies et al., 2014). Proper mitochondrial function, in turn, is essential to ensure the integrity of tendons and other connective tissues (Lowes et al., 2009; Thankam et al., 2018). With regard to rs268249346 (CHI 28: 18321523), the other genome-wide significant SNP for medial suspensory ligament, we found that it maps close to the ADAM metallopeptidase with thrombospondin type 1 motif 14 (ADAMTS14) locus 16 (Supplementary Table 3). This gene has been reported to encode a procollagen N- proteinase that cleaves the amino-propeptide of procollagen to allow the assembly of elongated and cylindrical collagen fibrils (Bekhouche and Colige, 2015), although one recent study indicates that the main role of this molecule would be the regulation of the immune response (Dupont et al., 2018). We also identified twelve SNPs showing significant associations with udder traits at the chromosome-wide level (Table 1, Supplementary Figures 4 and 7). The activating transcription factor 3 (ATF3) gene is located close to the SNP rs268273468 (CHI 16: 69617700) associated with medial suspensory ligament (Supplementary Table 3). Interestingly, this gene modulates the synthesis of collagen I and III (Zhou et al., 2011), and regulates matrix metalloproteinases which are fundamental for the development, renewal and remodeling of tendons (Guenzle et al., 2017). Another interesting association is that between the rs268288193 marker (CHI 19: 38362152) and teat placement. The region containing this SNP was associated with somatic cell count in Saanen goats (Martin et al., 2018). Interestingly, Lund et al., (1994) reported that cows with bad teat placement tend to be more susceptible to mastitis. Less than 1 Mb apart from this SNP, we identified the suppressor of cytokine signaling 7 (SOCS7) gene (Supplementary Table 3), which inhibits prolactin, growth hormone, and leptin signaling (Martens et al., 2005). Members of this gene family are regulators of mammary gland physiology. For instance, in mice SOCS1 and SOCS2 attenuate prolactin signaling, thus preventing premature lactation (Sutherland et al., 2007). Besides, SOCS3 is a key regulator of mammary gland involution (Sutherland et al., 2007). In dairy cattle, the polymorphism of the SOCS7 gene has been associated to different milk traits like protein yield and percentage, and milk yield (Arun et al., 2015). 17 The collagen type XIV α1 chain (COL14A1) gene is located just 0.1 Mb away from the rs268281312 SNP (CHI 14:868624) which is significantly associated to udder depth (Supplementary Table 3). This gene encodes a fibril-associated collagen that regulates fibrillogenesis ( Lindholm, 2019). In the human mammary gland, large amounts of type XIV collagen have been found in interlobular stroma, which, compared with intralobular stroma, contains densely packed collagen (Atherton et al., 1998). This differential distribution of type XIV collagen might have an important effect on the architecture of the mammary connective tissue. Identification of Genetic Determinants for Body Conformation Traits Although we did not detect any genome-wide significant association for body conformation traits, at the chromosome-wide level we found 19 SNPs significantly associated with angularity (4 SNPs), rump width (2 SNPs), rump angle (1 SNPs), chest width (5 SNPs), height (3 SNPs), body depth (3 SNPs) and bone quality (1 SNPs). These results are displayed in Table 2 as well as in Supplementary Figures 5 (Manhattan plots) and 8 (Q-Q plots). Several of the aforementioned SNPs map close to genes involved in bone homeostasis and skeletal development (Supplementary Table 3). For instance, the rs268245664 SNP which is associated with angularity maps to the parathyroid hormone 1 receptor (PTH1R) locus (Figure 5, Table 2, and Supplementary Table 3). This gene encodes a protein that acts as a receptor for parathyroid hormone (PHT) and parathyroid related peptide (PTHrP), two factors regulating mineral ion homeostasis (Mannstadt et al., 1999). The dysfunction of the PTH1R gene is associated with diseases that affect skeletal development and calcium homeostasis (Mannstadt et al., 1999). With regard to the association between angularity and the rs268265191 SNP (CHI 1: 35829812), it is 18 worth to mention that this marker maps near to the ephrin receptor A3 gene (EPHA3) (Supplementary Table 3). Ephrin receptors and their associated ligands are essential modulators of bone remodeling and they ensure an adequate coupling between bone resorption and formation (Edwards and Mundy, 2008). Interestingly, the inactivation of the ephrin-B1gene causes perinatal lethality, abdominal wall closure defects, and skeletal abnormalities, especially of the thoracic cage (Compagni et al., 2003). Another interesting gene is the one encoding the CGG triplet repeat binding protein 1 (CGGBP1), which is located 1 Mb away from rs268265191 (CHI 1: 35829812) (Supplementary Table 3). The polymorphism of this gene has been associated to several carcass traits in cattle (Calonge M., 2004), having a considerable effect on growth (Sevane et al., 2014). Moreover, Sevane et al. (2014) described a non- synonymous mutation, rs477676137 (c.206A>G, 1: 36060631), in the CGGBP1 gene associated with an increase in pelvis width and wither height measured at 9 months. On the other hand, the rs268255133 marker (CHI 18: 34112104), which is associated with rump width, maps close to the cadherin 11 gene (CDH11) (Supplementary Table 3),which modulates postnatal bone growth and osteoblast differentiation (Di Benedetto et al., 2010). Furthermore, the rs268262472 marker (CHI 27: 38084358), which is associated with rump angle, lies near to the spermatogenesis associated 4 gene (SPATA4), also involved in promoting osteoblast differentiation (Wang et al., 2011). Another marker of interest is rs268285858 (CHI 22: 45075519). This SNP was associated with body depth and it is located close to the Wnt family member 5A (WNT5A ) gene, (Supplementary Table 3), that regulates planar cell polarity signaling during embryonic development (Qian et al., 2007). The loss of the gene results in a shortened and widened cochlea in knockout mice embryos and severe shortening of the anterior-posterior axis and limb truncations due to abnormal convergent extension 19 (Yamaguchi et al., 1999; Qian et al., 2007; Andre et al., 2015). Finally, we would like to highlight that the rs268249930 (CHI 16: 26773653) marker associated with height co- localizes with the DNAH14 gene (Supplementary Table 3), which encodes an axonemal dynein heavy chain. Mutations in dynein genes can cause skeletal ciliopathies characterized by thoracic narrowing, short long bones and pelvis dysplasia (Yildiz, 2018). Identification of Genetic Determinants for Leg Structure Traits No genome-wide significant SNPs were found for leg structure traits, but 7 SNPs showed significant associations at the chromosome-wide level. These findings are reported in Table 3 as well as in Supplementary Figures 6 (Manhattan plot) and 9 (Q-Q plot). For the mobility trait, an association was found with the rs268236663 SNP (CHI 2:16211260), located 600 Kb away from the endothelin converting enzyme like 1 (ECEL1) gene (Supplementary Table 3), which encodes an endopeptidase member of the M13 family involved in the regulation of neuropeptide and peptide hormone activity. This molecule has an important function in the development of the neuromuscular junctions of the limbs in mice (Nagata et al., 2016). Mice lacking this gene display a poor arborization of the neuromotor nerves and a significant reduction of the number of neuromuscular junctions (Nagata et al., 2016). In humans, digital arthrogryposis is caused by mutations in ECEL1 and affected individuals show limited flexion of the knee and flexion of the fingers as well as muscular atrophy (Dieterich et al., 2013). Another interesting association was that between rear leg rear view and the genotype of the rs268286224 marker (CHI 24:42536581), which is positioned nearby the piezo type mechanosensitive ion channel component 2 (PIEZO2) gene, (Supplementary Table 3), 20 which is also involved in the etiology of digital arthrogryposis (Delle Vedove et al., 2016). Homozygous individuals for mutations inactivating PIEZO2 suffer from arthrogryposis and scoliosis (Haliloglu et al., 2017), while carriers of gain-of-function mutations can present multiple congenital contractures of limbs and variable absence of cruciate knee ligaments, amongst other symptoms (Coste et al., 2013). About the Genomic Architecture of Morphological Traits Comparing our results with those obtained by Martin et al., (2018) and by Mucha et al., (2018), we can state that there is a general lack of positional coincidence between the genomic regions associated with conformation traits in the aforementioned studies. Furthermore, in the study of Martin et al. (2018) different regions were identified in the Alpine and Saanen breeds as associated with body phenotypes, suggesting that this lack of concordance is not produced by technical factors. Indeed, we have found only 2 SNPs displaying genome-wide significant associations with morphological traits (2 SNPs). Such limited results are probably due to the fact that the inheritance of morphological traits in goats is highly polygenic, with many genetic variants with small effects determining phenotypic variation. In other words, the success of GWAS largely depends on the genomic architecture of the trait rather than on the magnitude of its heritability. Stature is a good example of this because it is a highly heritable and, at the same time, a highly polygenic trait. In cattle 163 genomic regions associated with stature have been detected, but they only explain 13.8% of the phenotypic variance (Bouwman et al., 2018). Similarly, in humans at least 180 genetic markers mostly segregating in populations of European descent explain ∼10% of the variance in height, a phenotype that has a large heritability close to ∼0.8 (Yang et al., 2010). This remarkable gap between genealogical heritability and the percentage of the phenotypic 21 variance explained by the SNPs is caused, at least in part, by the existence of hundreds or thousands of genetic determinants with small phenotypic effects on stature. A similar reasonment can be probably made for the majority of body, leg and udder morphological traits, because with the population sizes often employed in bovine, ovine or caprine GWAS generally the number of significant hits is very low (Schmid and Bennewitz, 2017). Currently, we do not know whether the lack of positional concordance between GWAS for morphology traits recorded in goats are due to the existence of a substantial genetic heterogeneity across populations or to the modest size of the populations employed in GWAS, a circumstance that limits the ability to detect variants with small effects that explain the majority of the phenotypic variance of morphological traits. Indeed, large GWAS analyses performed in humans have demonstrated that many of the genetic variants that are associated with height are shared between individuals of European and African descent (N’Diaye et al., 2011), thus reinforcing the notion that the magnitude of genetic heterogeneity between populations is greatly reduced when large sample sizes are used in GWAS. CONCLUSIONS The number of hits detected in our GWAS for body conformation, udder and leg traits was quite limited, a result that agrees well with those reported in previous studies (Martin et al., 2018; Mucha et al., 2018). Very probably, this outcome reflects the highly polygenic nature of morphological traits in ruminants. By comparing our results with previous reports, we have also detected a low positional concordance. This could be the consequence of genetic heterogeneity in the genetic determinism of 22 morphological traits or to the fact that GWAS carried out so far are too low powered to detect reliably the genetic determinants of such phenotypes. Despite these limitations, we have been able to identify several genes related with collagen synthesis (ATF3, ADAMTS14 and COL14A1), growth (CGGBP1), development (WNT5A and DNAH14), bone homeostasis and remodeling (PTH1R, CDH11, SPATA4, and EPHA3), limb development (ECEL1 and PIEZO2) and mammary physiology (SOCS7) mapping close to GWAS hits. Such information, combined with candidate gene sets generated in other GWAS, if possible, with much larger reference populations, could provide very valuable hints about the identity of the loci shaping the body, udder and leg morphology of goats. ACKNOWLEDGMENTS The authors are indebted to the Asociación Nacional de Criadores de Caprino de Raza Murciano-Granadina (Caprigran) for their cooperation in this study, and more specifically to Miguel García García and Teresa Novo Díaz, who collected all blood samples, and to Antonio Martín Ramírez who recorded all linear score measurements. This research was funded by the European Fund for Regional Development/Ministerio de Ciencia, Innovación y Universidades - Agencia Estatal de Investigación/ Project Reference (AGL2016-76108-R). We acknowledge the financial support from the Spanish Ministry of Economy and Competitiveness, through the “Severo Ochoa Programme for Centres of Excellence in R&D” 2016-2019 (SEV‐2015‐0533), and from the CERCA programme of the Generalitat de Catalunya. Emilio Mármol-Sánchez was funded with a FPU Ph.D. grant awarded by the Spanish Ministry of Education 23 (FPU15/01733). Maria Gracia Luigi-Sierra was funded with a FPI Ph.D. grant from the Spanish Ministry of Economy and Competitivity (BES-2017-079709). REFERENCES Andre, P., H. Song, W. Kim, A. Kispert, and Y. Yang. 2015. WNT5A and WNT11 regulate mammalian anterior-posterior axis elongation. Development 142:1516- 1527. doi:10.1242/dev.119065. Archer, S.C., M.J. Green, and J.N. Huxley. 2010. Association between milk yield and serial locomotion score assessments in UK dairy cows. J. Dairy Sci. 93:4045-4053. doi:10.3168/jds.2010-3062. Arun, S.J., P.C. Thomson, P.A. Sheehy, M.S. Khatkar, H.W. Raadsma, and P. Williamson. 2015. Targeted analysis reveals an important role of JAK-STAT-SOCS genes for milk production traits in Australian dairy cattle. Front. Genet. 6:342. doi:10.3389/fgene.2015.00342. Atherton, A.J., M.J. Warburton, M.J. O’Hare, P. Monaghan, D. Schuppan, and B.A. Gusterson. 1998. Differential expression of type XIV collagen/undulin by human mammary gland intralobular and interlobular fibroblasts. Cell Tissue Res. 291:507- 511. doi:10.1007/s004410051020. Aulchenko, Y.S., S. Ripke, A. Isaacs, and C.M. van Duijn. 2007. GenABEL: An R library for genome-wide association analysis. Bioinformatics. 23:1294-1296. doi:10.1093/bioinformatics/btm108. Bastin, C., and N. Gengler. 2013. Genetics of body condition score as an indicator of 24 dairy cattle fertility: A review. Biotechnol. Agron. Soc. Environ. 17:64-75. Bateman, A. 2019. UniProt: A worldwide hub of protein knowledge. Nucleic Acids Res. 47:D506–D515. doi:10.1093/nar/gky1049 Bekhouche, M., and A. Colige. 2015. The procollagen N-proteinases ADAMTS2, 3 and 14 in pathophysiology. Matrix Biol. 44-46:46-53. doi:10.1016/j.matbio.2015.04.001. Benjamini, Y. and Y. Hochberg. 1995. Controlling the false discovery rate - A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B. 57:289–300. doi:10.2307/2346101. Bickhart, D.M., B.D. Rosen, S. Koren, B.L. Sayre, A.R. Hastie, S. Chan, J. Lee, E.T. Lam, I. Liachko, S.T. Sullivan, J.N. Burton, H.J. Huson, J.C. Nystrom, C.M. Kelley, J.L. Hutchison, Y. Zhou, J. Sun, A. Crisà, F.A. Ponce De León, J.C. Schwartz, J.A. Hammond, G.C. Waldbieser, S.G. Schroeder, G.E. Liu, M.J. Dunham, J. Shendure, T.S. Sonstegard, A.M. Phillippy, C.P. Van Tassell, and T.P.L. Smith. 2017. Single-molecule sequencing and chromatin conformation capture enable de novo reference assembly of the domestic goat genome. Nat. Genet. 49:643–650. doi:10.1038/ng.3802. Boettcher, P.J., J.C.M. Dekkers, and B.W. Kolstad. 1998. Development of an udder health index for sire selection based on somatic cell score, udder conformation, and milking speed. J. Dairy Sci. 81:1157-1168. doi:10.3168/jds.S0022-0302(98)75678- 4. Bouwman, A.C., H.D. Daetwyler, A.J. Chamberlain, C.H. Ponce, M. Sargolzaei, F.S. Schenkel, G. Sahana, A. Govignon-Gion, S. Boitard, M. Dolezal, H. Pausch, R.F. 25 Brøndum, P.J. Bowman, B. Thomsen, B. Guldbrandtsen, M.S. Lund, B. Servin, D.J. Garrick, J. Reecy, J. Vilkki, A. Bagnato, M. Wang, J.L. Hoff, R.D. Schnabel, J.F. Taylor, A.A.E. Vinkhuyzen, F. Panitz, C. Bendixen, L.E. Holm, B. Gredler, C. Hozé, M. Boussaha, M.P. Sanchez, D. Rocha, A. Capitan, T. Tribout, A. Barbat, P. Croiseau, C. Drögemüller, V. Jagannathan, C. Vander Jagt, J.J. Crowley, A. Bieber, D.C. Purfield, D.P. Berry, R. Emmerling, K.U. Götz, M. Frischknecht, I. Russ, J. Sölkner, C.P. Van Tassell, R. Fries, P. Stothard, R.F. Veerkamp, D. Boichard, M.E. Goddard, and B.J. Hayes. 2018. Meta-analysis of genome-wide association studies for cattle stature identifies common genes that regulate body size in mammals. Nat. Genet. 50:362-367. doi:10.1038/s41588-018-0056-5. Calonge, M.E. 2004. Identificación de genes con expresión diferencial en tejido muscular de bovinos pertenecientes a los tres genotipos de la miostatina (mutación nt821(del11)). PhD Thesis. Universidad Complutense de Madrid, Spain. Campbell, J.R., and R.T. Marshall. 2016. Dairy production and processing: The science of milk and milk products. 1st ed. Waveland Press Inc., Long Grove, Illinois. Capote, J., A. Argüello, N. Castro, J.L. López, and G. Caja. 2006. Short communication: Correlations between udder morphology, milk yield, and milking ability with different milking frequencies in dairy goats. J. Dairy Sci. 89:2076– 2079. doi:10.3168/jds.S0022-0302(06)72276-7. Castañeda-Bustos, V.J., H.H. Montaldo, M. Valencia-Posadas, L. Shepard, S. Pérez- Elizalde, O. Hernández-Mendo, and G. Torres-Hernández. 2017. Linear and nonlinear genetic relationships between type traits and productive life in US dairy goats. J. Dairy Sci. 100:1232–1245. doi:10.3168/jds.2016-11313. 26 Chacón, E., F. Macedo, F. Velázquez, S.R. Paiva, E. Pineda, and C. McManus. 2011. Morphological measurements and body indices for Cuban Creole goats and their crossbreds. Rev. Bras. Zootec. 40:8. doi:10.1590/S1516-35982011000800007. Chang, C.C., C.C. Chow, L.C.A.M. Tellier, S. Vattikuti, S.M. Purcell, and J.J. Lee. 2015. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4:7. doi:10.1186/s13742-015-0047-8. Compagni, A., M. Logan, R. Klein, and R.H. Adams. 2003. Control of skeletal patterning by EphrinB1-EphB interactions. Dev. Cell. 5:217-230. doi:10.1016/S1534-5807(03)00198-9. Coste, B., G. Houge, M.F. Murray, N. Stitziel, M. Bandell, M.A. Giovanni, A. Philippakis, A. Hoischen, G. Riemer, U. Steen, V.M. Steen, J. Mathur, J. Cox, M. Lebo, H. Rehm, S.T. Weiss, J.N. Wood, R.L. Maas, S.R. Sunyaev, and A. Patapoutian. 2013. Gain-of-function mutations in the mechanically activated ion channel PIEZO2 cause a subtype of Distal Arthrogryposis. Proc. Natl. Acad. Sci. U.S.A. 110:4667-4672. doi:10.1073/pnas.1221400110. De la Fuente, L.F., C. Gonzalo, J.P. Sánchez, R. Rodríguez, J.A. Carriedo, and F.S. Primitivo. 2011. Genetic parameters of the linear body conformation traits and genetic correlations with udder traits, milk yield and composition, and somatic cell count in dairy ewes. Can. J. Anim. Sci. 91:585-591. doi:10.4141/cjas2010-031. Delgado, J.V., V. Landi, C.J. Barba, J. Fernández, M.M. Gómez, M.E. Camacho, M.A. Martínez, F.J. Navas, and J.M. León. 2018. Murciano-Granadina goat: A Spanish local breed ready for the challenges of the twenty-first century. Pages 205-219 in: Sustainable Goat Production in Adverse Environments: Volume II (J.Simões and 27 C. Gutiérrez, eds), Springer International Publishing, Heidelberg, Germany. Delle Vedove, A., M. Storbeck, R. Heller, I. Hölker, M. Hebbar, A. Shukla, O. Magnusson, S. Cirak, K.M. Girisha, M. O’Driscoll, B. Loeys, and B. Wirth. 2016. Biallelic loss of proprioception-related PIEZO2 causes muscular atrophy with perinatal respiratory distress, arthrogryposis, and scoliosis. Am. J. Hum. Genet. 99:1206-1216. doi:10.1016/j.ajhg.2016.09.019. Di Benedetto, A., M. Watkins, S. Grimston, V. Salazar, C. Donsante, G. Mbalaviele, G.L. Radice, and R. Civitelli. 2010. N-cadherin and cadherin 11 modulate postnatal bone growth and osteoblast differentiation by distinct mechanisms. J. Cell Sci. 123:2640-2648. doi:10.1242/jcs.067777. Dieterich, K., S. Quijano-Roy, N. Monnier, J. Zhou, J. Fauré, D.A. Smirnow, R. Carlier, C. Laroche, P. Marcorelles, S. Mercier, A. Mégarbané, S. Odent, N. Romero, D. Sternberg, I. Marty, B. Estournet, P.S. Jouk, J. Melki, and J. Lunardi. 2013. The neuronal endopeptidase ECEL1 is associated with a distinct form of recessive distal arthrogryposis. Hum. Mol. Genet. 22:1483-1492. doi:10.1093/hmg/dds514. Dupont, L., G. Ehx, M. Chantry, C. Monseur, C. Leduc, L. Janssen, D. Cataldo, M. Thiry, C. Jerome, J.M. Thomassin, B. Nusgens, J. Dubail, F. Baron, and A. Colige. 2018. Spontaneous atopic dermatitis due to immune dysregulation in mice lacking ADAMTS2 and 14. Matrix Biol. 70:140-157. doi:10.1016/j.matbio.2018.04.002. Edwards, C.M., and G.R. Mundy. 2008. Eph receptors and ephrin signaling pathways: A role in bone homeostasis. Int. J. Med. Sci. 5:263-272. doi:10.7150/ijms.5.263. Egger-Danner, C., J.B. Cole, J.E. Pryce, N. Gengler, B. Heringstad, A. Bradley, and K.F. Stock. 2014. Invited review: Overview of new traits and phenotyping 28 strategies in dairy cattle with a focus on functional traits. Animal. 9:191–207. doi:10.1017/S1751731114002614. Guan, D., V. Landi, M.G. Luigi-Sierra, J.V. Delgado, X. Such, A. Castelló, B. Cabrera, E. Mármol-Sánchez, J. Fernández-Alvarez, J.L.R. de la Torre Casañas, A. Martínez, J. Jordana, and M. Amills. 2020. Analyzing the genomic and transcriptomic architecture of milk traits in Murciano-Granadina goats. J. Anim. Sci. Biotechnol. 11:35. doi:10.1186/s40104-020-00435-4. Guarín, J.F., M.G. Paixão, and P.L. Ruegg. 2017. Association of anatomical characteristics of teats with quarter-level somatic cell count. J. Dairy Sci. 100:643- 652. doi:10.3168/jds.2016-11459. Guenzle, J., L.J. Wolf, N.W.C. Garrelfs, J.M. Goeldner, N. Osterberg, C.R. Schindler, J.E. Saavedra, and A. Weyerbrock. 2017. ATF3 reduces migration capacity by regulation of matrix metalloproteinases via NFκB and STAT3 inhibition in glioblastoma. Cell Death Discov. 3:17006. doi:10.1038/cddiscovery.2017.6. Haliloglu, G., K. Becker, C. Temucin, B. Talim, N. Küçüksąhin, M. Pergande, S. Motameny, P. Nürnberg, U. Aydingoz, H. Topaloglu, and S. Cirak. 2017. Recessive PIEZO2 stop mutation causes distal arthrogryposis with distal muscle weakness, scoliosis and proprioception defects. J. Hum. Genet. 62:497-501. doi:10.1038/jhg.2016.153. Huang, D.W., B.T. Sherman, and R.A. Lempicki. 2009. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4:44–57. doi:10.1038/nprot.2008.211. Kinsella, R.J., A. Kähäri, S. Haider, J. Zamora, G. Proctor, G. Spudich, J. Almeida- 29 King, D. Staines, P. Derwent, A. Kerhornou, P. Kersey, and P. Flicek. 2011. Ensembl BioMarts: A hub for data retrieval across taxonomic space. Database. 2011. doi:10.1093/database/bar030. Legarra, A., and E. Ugarte. 2005. Genetic parameters of udder traits, somatic cell score, and milk yield in Latxa sheep. J. Dairy Sci. 88:2238-2245. doi:10.3168/jds.S0022- 0302(05)72899-X. Lindholm M., T. Manon-Jensen and M.A. Karsdal. 2019. Type XIV collagen. Pages 121-125 in Biochemistry of Collagens, Laminins and Elastin. Vol. 2. (M.A. Karsdal, ed.) Academic Press, London, UK. Lowes, D.A., C. Wallace, M.P. Murphy, N.R. Webster, and H.F. Galley. 2009. The mitochondria targeted antioxidant MitoQ protects against fluoroquinolone-induced oxidative stress and mitochondrial membrane damage in human Achilles tendon cells. Free Radic. Res. 43:323-328. doi:10.1080/10715760902736275. Lund, T., F. Miglior, J.C.M. Dekkers, and E.B. Burnside. 1994. Genetic relationships between clinical mastitis, somatic cell count, and udder conformation in Danish Holsteins. Livest. Prod. Sci. 3:243-251. doi:10.1016/0301-6226(94)90203-8. Manfredi, E., A. Piacere, P. Lahaye, and V. Ducrocq. 2001. Genetic parameters of type appraisal in Saanen and Alpine goats. Livest. Prod. Sci. 70:183–189. doi:10.1016/S0301-6226(01)00180-4. Mannstadt, M., H. Jüppner, and T.J. Gardella. 1999. Receptors for PTH and PTHrP: Their biological importance and functional properties. Am. J. Physiol. 277:F665- F675. doi: 10.1152/ajprenal.1999.277.5.F665. 30 Martens, N., G. Uzan, M. Wery, R. Hooghe, E.L. Hooghe-Peters, and A. Gertler. 2005. Suppressor of cytokine signaling 7 inhibits prolactin, growth hormone, and leptin signaling by interacting with STAT5 or STAT3 and attenuating their nuclear translocation. J. Biol. Chem. 280:13817-13823. doi:10.1074/jbc.M411596200. Martin, P., I. Palhière, C. Maroteau, V. Clément, I. David, G.T. Klopp, and R. Rupp. 2018. Genome-wide association mapping for type and mammary health traits in French dairy goats identifies a pleiotropic region on chromosome 19 in the Saanen breed. J. Dairy Sci. 101:5214–5226. doi:10.3168/jds.2017-13625. McLaren, A., S. Mucha, R. Mrode, M. Coffey, and J. Conington. 2016. Genetic parameters of linear conformation type traits and their relationship with milk yield throughout lactation in mixed-breed dairy goats. J. Dairy Sci. 99:5516–5525. doi:10.3168/jds.2015-10269. Miglior, F., A. Fleming, F. Malchiodi, L.F. Brito, P. Martin, and C.F. Baes. 2017. A 100-Year Review: Identification and genetic selection of economically important traits in dairy cattle. J. Dairy Sci. 100:10251–10271. doi:10.3168/jds.2017-12968. Miller, S.A., D.D. Dykes, and H.F. Polesky. 1988. A simple salting out procedure for extracting DNA from human nucleated cells. Nucleic Acids Res. 16:1215. doi:10.1093/nar/16.3.1215. Montaldo, H.H., and E. Manfredi. 2002. Organisation of selection programmes for dairy goats. No. 01–35 in 7th World Congr. Genet. Appl. Livest. Prod. Montpellier, France. Mucha, S., R. Mrode, M. Coffey, M. Kizilaslan, S. Desire, and J. Conington. 2018. Genome-wide association study of conformation and milk yield in mixed-breed 31 dairy goats. J. Dairy Sci. 101:2213–2225. doi:10.3168/jds.2017-12919. N’Diaye, A., G.K. Chen, C.D. Palmer, B. Ge, B. Tayo, R.A. Mathias, J. Ding, M.A. Nalls, A. Adeyemo, V. Adoue, C.B. Ambrosone, L. Atwood, E. V. Bandera, L.C. Becker, S.I. Berndt, L. Bernstein, W.J. Blot, E. Boerwinkle, A. Britton, G. Casey, S.J. Chanock, E. Demerath, S.L. Deming, W.R. Diver, C. Fox, T.B. Harris, D.G. Hernandez, J.J. Hu, S.A. Ingles, E.M. John, C. Johnson, B. Keating, R.A. Kittles, L.N. Kolonel, S.B. Kritchevsky, L. Marchand, K. Lohman, J. Liu, R.C. Millikan, A. Murphy, S. Musani, C. Neslund-Dudas, K.E. North, S. Nyante, A. Ogunniyi, E.A. Ostrander, G. Papanicolaou, S. Patel, C.A. Pettaway, M.F. Press, S. Redline, J.L. Rodriguez-Gil, C. Rotimi, B.A. Rybicki, B. Salako, P.J. Schreiner, L.B. Signorello, A.B. Singleton, J.L. Stanford, A.H. Stram, D.O. Stram, S.S. Strom, B. Suktitipat, M.J. Thun, J.S. Witte, L.R. Yanek, R.G. Ziegler, W. Zheng, X. Zhu, J.M. Zmuda, A.B. Zonderman, M.K. Evans, Y. Liu, D.M. Becker, R.S. Cooper, T. Pastinen, B.E. Henderson, J.N. Hirschhorn, G. Lettre, and C.A. Haiman. 2011. Identification, replication, and fine-mapping of loci associated with adult height in individuals of African ancestry. PLoS Genet. 7: e1002298. doi:10.1371/journal.pgen.1002298. Nagata, K., S. Kiryu-Seo, H. Tamada, F. Okuyama-Uchimura, H. Kiyama, and T.C. Saido. 2016. ECEL1 mutation implicates impaired axonal arborization of motor nerves in the pathogenesis of distal arthrogryposis. Acta Neuropathol. 132:111- 126. doi:10.1007/s00401-016-1554-0. Paradies, G., V. Paradies, V. De Benedictis, F.M. Ruggiero, and G. Petrosillo. 2014. Functional role of cardiolipin in mitochondrial bioenergetics. Biochim. Biophys. Acta. 1837:408-417. doi:10.1016/j.bbabio.2013.10.006. 32 Pérez-Cabal, M.A., E. Legaz, I. Cervantes, L.F. de la Fuente, R. Martínez, F. Goyache, and J.P. Gutiérrez. 2013. Association between body and udder morphological traits and dairy performance in Spanish Assaf sheep. Arch. Anim. Breed. 56:430-442. doi:10.7482/0003-9438-56-042. Pugh, D.G., and A.N. Baird. 2012. Sheep and Goat Medicine. 2nd ed. Elsevier, St. Louis, Missouri, USA. Qian, D., C. Jones, A. Rzadzinska, S. Mark, X. Zhang, K.P. Steel, X. Dai, and P. Chen. 2007. WNT5A functions in planar cell polarity regulation in mice. Dev. Biol. 306:121-133. doi:10.1016/j.ydbio.2007.03.011. R Core team. 2017. R: A Language and environment for statistical computing. R Found. Stat. Comput. Vienna, Austria. ISBN 3-900051-07-0, URL http//www.R- project.org/. Rogers, G.W., G.L. Hargrove, T.J. Lawlor, and J.L. Ebersole. 1991. Correlations among linear type traits and somatic cell counts. J. Dairy Sci. 74:1087-1091. doi:10.3168/jds.S0022-0302(91)78259-3. Rupp, R., and D. Boichard. 1999. Genetic parameters for clinical mastitis, somatic cell score, production, udder type traits, and milking ease in first lactation Holsteins. J. Dairy Sci. 82:2198-2204. doi:10.3168/jds.S0022-0302(99)75465-2. Rupp, R., and D. Boichard. 2003. Genetics of resistance to mastitis in dairy cattle. Vet. Res. 34:671–688. doi:10.1051/vetres:2003020. Rupp, R., V. Clément, A. Piacere, C. Robert-Granié, and E. Manfredi. 2011. Genetic parameters for milk somatic cell score and relationship with production and udder 33 type traits in dairy Alpine and Saanen primiparous goats. J. Dairy Sci. 7: 3629- 3634. doi:10.3168/jds.2010-3694. Sánchez-Rodríguez, M. 2012. Valoración morfológica del ganado caprino lechero. Juzgamiento y calificación. 2nd Ed. Editorial Servet, Zaragoza, Spain. Schmid, M., and J. Bennewitz. 2017. Invited review: Genome-wide association analysis for quantitative traits in livestock - A selective review of statistical models and experimental designs. Arch. Anim. Breed. 94:3629-3634. doi:10.5194/aab-60-335- 2017. Sevane, N., E. Armstrong, P. Wiener, R. Pong Wong, S. Dunner, V. Amarger, D. Delourme, H. Levéziel, S. Boitard, B. Mangin, J. Cañón, M.L. Checa, D. García, M.E. Miranda, R. Pérez, M. Christensen, P. Ertbjerg, A. Crisá, C. Marchitelli, A. Valentini, S. Failla, S. Gigli, J.F. Hocquette, G. Nute, I. Richardson, J.L. Olleta, B. Panea, C. Sañudo, N. Razzaq, G. Renand, and J.L. Williams. 2014. Polymorphisms in twelve candidate genes are associated with growth, muscle lipid profile and meat quality traits in eleven European cattle breeds. Mol. Biol. Rep. 41:4721– 4731. doi:10.1007/s11033-014-3343-y. Seykora, A.J., and B.T. McDaniel. 1985. Udder and teat morphology related to mastitis resistance: A review. J. Dairy Sci. 68:2087–2093. doi:10.3168/jds.S0022- 0302(85)81072-9. Shelton, M. 1978. Reproduction and breeding of goats. J. Dairy Sci. 61:994–1010. doi:10.3168/jds.S0022-0302(78)83680-7. Shim, H., D.I. Chasman, J.D. Smith, S. Mora, P.M. Ridker, D.A. Nickerson, R.M. Krauss, and M. Stephens. 2015. A multivariate genome-wide association analysis 34 of 10 LDL subfractions, and their response to statin treatment, in 1868 Caucasians. PLoS One 10:e0120758. doi:10.1371/journal.pone.0120758. Sutherland, K.D., G.J. Lindeman, and J.E. Visvader. 2007. Knocking off SOCS genes in the mammary gland. Cell Cycle 6:799-803. doi:10.4161/cc.6.7.4037. Thankam, F.G., I.S. Chandra, A.N. Kovilam, C.G. Diaz, B.T. Volberding, M.F. Dilisio, M.M. Radwan, R.M. Gross, and D.K. Agrawal. 2018. Amplification of mitochondrial activity in the healing response following rotator cuff tendon injury. Sci. Rep. 8:17027. doi:10.1038/s41598-018-35391-7. Wang, X., K. Harimoto, J. Liu, J. Guo, S. Hinshaw, Z. Chang, and Z. Wang. 2011. SPATA4 promotes osteoblast differentiation through Erk-activated Runx2 pathway. J. Bone Miner. Res. 8:1964-1973. doi:10.1002/jbmr.394. Yamaguchi, T.P., A. Bradley, A.P. McMahon, and S. Jones. 1999. A WNT5A pathway underlies outgrowth of multiple structures in the vertebrate embryo. Development 126:1211-1223. Yang, J., B. Benyamin, B.P. McEvoy, S. Gordon, A.K. Henders, D.R. Nyholt, P.A. Madden, A.C. Heath, N.G. Martin, G.W. Montgomery, M.E. Goddard, and P.M. Visscher. 2010. Common SNPs explain a large proportion of the heritability for human height. Nat. Genet. 42:565-569. doi:10.1038/ng.608. Yang, Y., J. Cao, and Y. Shi. 2004. Identification and characterization of a gene encoding human LPGAT1, an endoplasmic reticulum-associated lysophosphatidylglycerol acyltransferase. J. Biol. Chem. 279:55866-55874. doi:10.1074/jbc.M406710200. 35 Yildiz, A. 2018. Single-molecule dynein motor mechanics in vitro. Pages 113-135 in Dynein Mechanics, Dysfunction, and Disease. 2nd ed. S.M. King, ed. Academic Press, London, UK. Zhou, H., D.F. Shen, Z.Y. Bian, J. Zong, W. Deng, Y. Zhang, Y.Y. Guo, H. Li, and Q.Z. Tang. 2011. Activating transcription factor 3 deficiency promotes cardiac hypertrophy, dysfunction, and fibrosis induced by pressure overload. PLoS One. 6: e26744. doi:10.1371/journal.pone.0026744. Zhou, X., and M. Stephens. 2012. Genome-wide efficient mixed-model analysis for association studies. Nat. Genet. 44:821–824. doi:10.1038/ng.2310. Zujovic, M., N. Memisi, V. Bogdanovic, and Z. Tomic. 2011. Correlation between body measurements and milk production of goats in different lactations. Biotechnol. Anim. Husb. 27:217-225. doi:10.2298/bah1102217z. 36 Table 1. Genome-wide and chromosome-wide significant associations between SNPs and udder traits (MSL, medial suspensory ligament; TD, teat diameter; TP, teat placement; UD, udder depth; UW, udder width) recorded in 811 Murciano-Granadina goats. Trait Chr1 rs2 Pos3 A14 MAF5 β ± SE6 P-value7 q-value8 PVE9 16 rs268273468 69617700 C 0.469 -0.220 ± 0.045 1.45E-06 0.034 0.010 Genome Wide MSL 28 rs268249346 18321523 A 0.485 -0.241 ± 0.049 1.28E-06 0.034 0.008 MSL 16 rs268273468 69617700 C 0.469 -0.220 ± 0.045 1.45E-06 0.002 0.010 28 rs268249346 18321523 A 0.485 -0.241 ± 0.049 1.28E-06 0.001 0.008 TD 28 rs268248647 9579785 G 0.460 -0.214 ± 0.054 6.93E-05 0.042 0.007 Chromosome Wide 28 rs268243765 10586645 G 0.404 -0.215 ± 0.055 9.45E-05 0.042 0.007 TP 9 rs268282545 20855427 G 0.498 -0.208 ± 0.047 1.09E-05 0.020 0.009 19 rs268288193 38362152 G 0.181 0.243 ± 0.059 4.05E-05 0.047 0.006 25 rs268246864 40499453 G 0.391 -0.205 ± 0.047 1.32E-05 0.011 0.009 26 rs268291440 37970238 A 0.316 -0.238 ± 0.052 5.33E-06 0.005 0.008 UD 14 rs268281312 868624 A 0.263 0.216 ± 0.051 2.83E-05 0.041 0.008 14 rs268276674 90209603 A 0.480 0.193 ± 0.047 4.48E-05 0.041 0.009 15 rs268268821 47457448 A 0.430 -0.207 ± 0.050 3.14E-05 0.049 0.008 UW 27 rs268251218 10045474 G 0.444 0.277 ± 0.062 1.00E-05 0.009 0.005 1 Chr, chromosome; 2rs, identifier code of the SNP according to the RefSNP database; 3 Pos, position in base pairs; 4A1, minority allele; 5MAF, allele frequency; 6β ± SE, allelic substitution effect ± standard error; 7P-value, raw P-values; 8q-value, P-values corrected for multiple testing using a false discovery rate approach; 9PVE, percentage of proportion of variance in phenotype explained by a given SNP. 37 Table 2. Chromosome-wide significant associations between SNPs and body traits (ANG, angularity; BD, body depth; BQ, bone quality; CW, chest width; HT, height; RA, rump angle; RW, rump width) recorded in 811 Murciano-Granadina goats. Trait Chr1 rs2 Pos3 A14 MAF5 β ± SE6 P-value7 q- PVE9 value8 ANG 1 rs268265191 35829812 G 0.496 -0.199 ± 0.045 0.032 4.980 0.010 1 rs268280713 135255901 A 0.352 -0.200 ± 0.048 0.049 4.504 0.009 22 rs268245662 52518237 A 0.074 -0.358 ± 0.083 0.010 4.759 0.003 22 rs268245664 52649718 G 0.043 -0.512 ± 0.105 0.001 5.900 0.002 BD 9 rs268236956 24849991 G 0.281 0.224 ± 0.050 0.012 5.161 0.008 13 rs268290607 55121535 A 0.394 -0.190 ± 0.045 0.044 4.553 0.010 Chromosome Wide 22 rs268285858 45075519 G 0.48 -0.180 ± 0.044 0.045 4.393 0.011 BQ 15 rs268289470 13689629 A 0.394 0.217 ± 0.052 0.012 5.109 0.008 CW 3 rs268243320 53360997 C 0.139 -0.265 ± 0.062 0.044 4.715 0.005 7 rs268240258 51615185 A 0.311 -0.214 ± 0.045 0.006 5.548 0.010 17 rs268248977 26865599 A 0.454 0.174 ± 0.040 0.021 4.824 0.013 17 rs268248975 26942955 G 0.252 0.193 ± 0.049 0.044 4.022 0.008 17 rs268264934 58534119 G 0.494 0.171 ± 0.041 0.027 4.411 0.012 HT 16 rs268249920 26376843 A 0.419 0.172 ± 0.043 0.031 4.213 0.011 16 rs268249930 26773653 G 0.348 -0.207 ± 0.044 0.005 5.477 0.011 16 rs268236696 59779757 A 0.126 -0.238 ± 0.059 0.031 4.262 0.006 RA 27 rs268262472 38084358 G 0.207 -0.255 ± 0.063 0.044 4.299 0.005 RW 17 rs268264930 58698058 G 0.359 0.182 ± 0.040 0.009 5.190 0.013 18 rs268255133 34112104 G 0.134 -0.259 ± 0.055 0.003 5.553 0.007 1 Chr, chromosome; 2rs, identifier code of the SNP according to the RefSNP database; 3 Pos, position in base pairs; 4A1, minority allele; 5MAF, allele frequency; 6β ± SE, allelic substitution effect ± standard error; 7P-value, raw P-values; 8q-value, P-values corrected for multiple testing using a false discovery rate approach; 9PVE, percentage of proportion of variance in phenotype explained by a given SNP. 38 Table 3. Chromosome-wide significant associations between SNPs and leg traits (MOB, mobility; RLR, rear legs rear view; RLS, rear legs side view) recorded in 811 Murciano-Granadina goats. Trait Chr1 rs2 Pos3 A14 MAF5 β ± SE6 P-value7 q- PVE9 value8 MOB 2 rs268236663 16211260 G 0.237 0.325 ± 8.13E-06 0.022 0.004 0.072 RLR 17 rs268258221 21284278 G 0.367 -0.255 2.98E-05 0.041 0.006 ± 0.061 Chromosome Wide 19 rs268288174 39252784 G 0.188 0.340 ± 4.69E-06 0.005 0.004 0.074 24 rs268259571 38761188 A 0.422 -0.278 1.78E-05 0.023 0.005 ± 0.065 24 rs268286224 42536581 A 0.117 -0.376 4.26E-05 0.027 0.002 ± 0.092 RLS 3 rs268254620 8818355 G 0.356 -0.188 4.41E-05 0.050 0.010 ± 0.046 3 rs268285963 43496497 A 0.3 0.211 ± 1.20E-05 0.027 0.009 0.048 1 Chr, chromosome; 2rs, identifier code of the SNP according to the RefSNP database; 3 Pos, position in base pairs; 4A1, minority allele; 5MAF, allele frequency; 6β ± SE, allelic substitution effect ± standard error; 7P-value, raw P-values; 8q-value, P-values corrected for multiple testing using a false discovery rate approach; 9PVE, proportion of variance in phenotype explained by a given SNP. 39 LEGENDS TO FIGURES Figure 1A. Negative log10 P-values (y-axis) of the associations between SNPs and the medial suspensory ligament phenotype are plotted against the genomic location of each SNP marker (x-axis). Markers on different chromosomes are denoted by different colors. The blue line indicates the significance of the association after false discovery rate correction for multiple testing (q-value= 0.05). 1B. Quantile-Quantile plots corresponding to the genome-wide corrected P-values of the GWAS for the trait medial suspensory ligament and its lambda (λ) inflation factor. Figure 2A. Heatmap depicting Pearson correlations between morphological and milk yield phenotypes recorded in 811 Murciano-Granadina goats. The intensity of the color indicates the magnitude of the correlation 2B. Heatmap displaying the significances of the Pearson correlations shown in 2A. The intensity of the color indicates the significance of the association. The following abbreviations have been used, HT, height; CW, chest width; BD, body depth; RW, rump width; RA, rump angle; ANG, angularity; BQ, bone quality; RLR, rear legs rear view; RLS, rear legs side view; MOB, mobility; FUA, fore udder attachment; RUH, rear udder height; MSL, medial suspensory ligament; UW, udder width; UD, udder depth; TP, teat placement; TD, teat diameter; RMY, total milk yield; SCC, somatic cell count; MY210, milk yield corrected at 210 days milking; MY240, milk yield corrected at 240 days milking; M305, milk yield corrected at 305 days milking. 40 41 42 Supplementary Table 1. Descriptive statistics of body, udder and leg phenotypes recorded in 811 Murciano-Granadina goats and expressed as linear scores. Phenotype Mean ± SE1 Min2 Max3 Height 3.205 ± 0.092 1 8 Chest width 4.857 ± 0.087 2 9 Body depth 5.821 ± 0.061 3 9 Rump width 5.167 ± 0.094 3 8 Rump angle 5.655 ± 0.090 4 7 Angularity 5.204 ± 0.086 2 9 Bone quality 7.421 ± 0.094 5 9 Fore udder 5.917 ± 0.089 2 8 attachment Rear udder height 6.365 ± 0.084 3 8 Medial suspensory 3.212 ± 0.090 1 8 ligament Udder width 7.147 ± 0.098 3 9 Udder depth 4.406 ± 0.095 1 8 Teat placement 6.536 ± 0.085 3 9 Teat diameter 5.152 ± 0.094 2 9 Rear legs side view 6.673 ± 0.079 4 8 Rear legs rear view 3.342 ± 0.083 2 6 Mobility 6.926 ± 0.075 4 8 1 Mean ± SE, mean ± standard error;2Min, minimum values measured in the population; 3 Max, maximum values measured in the population. 43 1 Supplementary Table 2. Phenotypic Pearson correlations matrix among morphological and milk yield traits recorded in 811 Murciano- 2 Granadina goats. The lower diagonal shows the correlation between each pair of traits, the upper diagonal indicates the P-value of each 3 correlation. 4 Attached in excel file. 5 6 7 8 9 10 11 1 12 Supplementary Table 3. Distance in base pairs between SNPs showing significant 13 associations with morphological traits and the closest functional and positional 14 candidate gene identified in the goat ARS1 reference genome. Trait Gene Chr1 Pos1 2 Pos23 rs4 Pos5 Distance Gene- SNP6 MSL7 LPGAT1 16 70237593 70372062 rs268273468 69617700 619893 MSL7 ATF3 16 69579015 69584224 rs268273468 69617700 38685 MSL7 ADAMTS14 28 18733706 18830075 rs268249346 18321523 412183 TP8 SOCS7 19 38654621 38680473 rs268288193 38362152 292469 UD9 COL14A1 14 490458 701160 rs268281312 868624 167464 ANG10 PTH1R 22 52630037 52650450 rs268245664 52649718 19681 RW11 CDH11 18 33975536 34062591 rs268255133 34112104 136568 RA12 SPATA4 27 37814782 37824810 rs268262472 38084358 269576 ANG10 EPHA3 1 35895723 36304923 rs268265191 35829812 65911 ANG10 CGGBP1 1 34861253 34861756 rs268265191 35829812 968559 BD13 WNT5A 22 45489607 45511154 rs268285858 45075519 414088 HT14 DNAH14 16 26597510 26999792 rs268249930 26773653 - MOB15 ECEL1 2 15610337 15619144 rs268236663 16211260 592116 RLR16 PIEZO2 24 42556947 42810419 rs268286224 42536581 20366 1 15 Chr, chromosome; 2Pos1, start position of the gene in base pairs; 3Pos2, end position 16 of the gene in base pairs; 4rs, identifier code of the SNP according to the RefSNP 17 database; 5Pos, position of the SNP; 6Distance Gene-SNP, distance in base pairs 18 between the significant related SNP and the functional putative gene; 7MSL, medial 19 suspensory; 8TP, teat placement; 9UD, udder depth; 10ANG, angularity; 11RW, rump 20 width; 12RA, rump angle; 13BD, body depth; 14HT, height; 15MOB, mobility; 16RLR, 21 rear legs rear view. 22 . 23 24 25 1 26 Supplementary Figure 1. Principal component analysis based on the Illumina Goat 27 SNP50 BeadChip (Illumina inc, San Diego, CA) genotypes of 811 Murciano-Granadina 28 distributed in 13 farms (each one encoded by a different color) and used in the current 29 GWAS. 30 31 32 33 34 35 36 37 38 39 40 2 41 Supplementary Figure 2. Histograms depicting the number of Murciano-Granadina 42 goats corresponding to each one of the linear score categories defining udder traits. The 43 following abbreviations have been used, FUA, fore udder attachment; RUH, rear udder 44 height; MSL, medial suspensory ligament; UW, udder width; UD, udder depth; TP, teat 45 placement; TD, teat diameter. 46 47 48 49 3 50 Supplementary Figure 3. Histograms depicting the number of Murciano-Granadina 51 goats corresponding to each one of the linear score categories defining body (A) and leg 52 (B) traits. The following abbreviations have been used, HT, height; CW, chest width; 53 BD, body depth; RW, rump width; RA, rump angle; ANG, angularity; BQ, bone 54 quality; RLR, rear legs rear view; RLS, rear legs side view; MOB, mobility. 55 4 56 Supplementary Figure 4. Manhattan plots depicting negative log10 P-values (y-axis) of 57 the associations at the chromosome-wide level between SNPs and udder morphology 58 traits plotted against the genomic location of each SNP marker (x-axis). Only 59 chromosomes hosting significant associations are shown (q-value < 0.05). The 60 following abbreviations have been used, MSL, medial suspensory ligament; TD, teat 61 diameter; TP, teat placement; UD, udder depth; UW, udder width. 62 63 5 64 Supplementary Figure 5. Manhattan plots depicting negative log10 P-values (y-axis) of 65 the associations at the chromosome-wide level between SNPs and body morphology 66 traits plotted against the genomic location of each SNP marker (x-axis). Only 67 chromosome hosting significant associations are shown (q-value < 0.05).The following 68 abbreviations have been used, ANG, angularity; BD, body depth; BQ, bone quality; 69 CW, chest width; HT, height; RA, rump angle; RW, rump width. 70 6 71 Supplementary Figure 6. Manhattan plots depicting negative log10 P-values (y-axis) of 72 the associations at the chromosome-wide level between SNPs and leg traits plotted 73 against the genomic location of each SNP marker (x-axis). Only chromosomes hosting 74 significant associations are shown (q-value < 0.05). The following abbreviations have 75 been used, MOB, mobility; RLR, rear legs rear view; RLS, rear legs side view. 76 77 78 79 80 81 82 83 84 85 86 87 7 88 Supplementary Figure 7. Quantile-Quantile plots showing the expected distribution of 89 the -log10 P-values (x-axis) compared to the observed -log10 P-values (y-axis) in the 90 GWAS performed at the chromosome wide level for udder traits.The mean lambda (λ) 91 inflation factor was 1.060 ± 0.201. The following abbreviations have been used, MSL, 92 medial suspensory ligament; TD, teat diameter; TP, teat placement; UD, udder depth; 93 UW, udder width. 94 8 95 Supplementary Figure 8. Quantile-Quantile plots showing the expected distribution of 96 the -log10 P-values (x-axis) compared to the observed -log10 P-values (y-axis) in the 97 GWAS performed at the chromosome wide level for body traits. The mean lambda (λ) 98 inflation factor was 1.090 ± 0.113. The following abbreviations have been used, ANG, 99 angularity; BD, body depth; BQ, bone quality; CW, chest width; HT, height; RA, rump 100 angle; RW, rump width. 101 9 102 Supplementary Figure 9. Quantile-Quantile plots showing the expected distribution of 103 the -log10 P-values (x-axis) compared to the observed -log10 P-values (y-axis) in the 104 GWAS performed at the chromosome wide level for leg traits goats. The mean lambda 105 (λ) inflation factor was 1.098 ± 0.272. The following abbreviations have been used, 106 MOB, mobility; RLR, rear legs rear view; RLS, rear legs side view. 107 108 109 110 111 112 113 114 115 116 117 10