Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 https://doi.org/10.1186/s40104-022-00684-5 RESEARCH Open Access Genomic patterns of homozygosity and inbreeding depression in Murciano- Granadina goats María Gracia Luigi-Sierra1, Almudena Fernández2, Amparo Martínez3, Dailu Guan1, Juan Vicente Delgado3, Javier Fernández Álvarez3, Vincenzo Landi4, Francesc Xavier Such5, Jordi Jordana6, María Saura2* and Marcel Amills1,6* Abstract Background: Inbreeding depression can adversely affect traits related to fitness, reproduction and productive performance. Although current research suggests that inbreeding levels are generally low in most goat breeds, the impact of inbreeding depression on phenotypes of economic interest has only been investigated in a few studies based on genealogical data. Results: We genotyped 1040 goats with the Goat SNP50 BeadChip. This information was used to estimate different molecular inbreeding coefficients and characterise runs of homozygosity and homozygosity patterns. We detected 38 genomic regions with increased homozygosity as well as 8 ROH hotspots mapping to chromosomes 1, 2, 4, 6, 14, 16 and 17. Eight hundred seventeen goats with available records for dairy traits were analysed to evaluate the potential consequences of inbreeding depression on milk phenotypes. Four regions on chromosomes 8 and 25 were significantly associated with inbreeding depression for the natural logarithm of the somatic cell count. Notably, these regions contain several genes related with immunity, such as SYK, IL27, CCL19 and CCL21. Moreover, one region on chromosome 2 was significantly associated with inbreeding depression for milk yield. Conclusions: Although genomic inbreeding levels are low in Murciano-Granadina goats, significant evidence of inbreeding depression for the logarithm of the somatic cell count, a phenotype closely associated with udder health and milk yield, have been detected in this population. Minimising inbreeding would be expected to augment economic gain by increasing milk yield and reducing the incidence of mastitis, which is one of the main causes of dairy goat culling. Keywords: Goat, Inbreeding, Milk yield, Murciano-Granadina, Somatic cell score * Correspondence: saura.maria@inia.es; marcel.amills@uab.cat 2 Departamento de Mejora Genética Animal, INIA, Carretera de la Coruña km 7,5, 28040 Madrid, Spain 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain Full list of author information is available at the end of the article © The Author(s). 2022 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 2 of 14 Background Several studies have used genomic methods to deter- Inbreeding is defined as the mating of individuals that mine the levels of inbreeding in goat populations with a are related to each other more closely than the average broad geographic distribution [19–21]. A recent investi- relationship within the concerned population [1]. In live- gation carried out by Bertolini et al. [19] revealed that stock species, the magnitude of inbreeding has been short ROH (< 3 Mb) are particularly abundant in world- traditionally measured through genealogical information wide goat populations. Moreover, five regions on caprine [2]. However, pedigree-based estimates are affected by chromosomes (CHI) 11, 12, and 18 contained ROH hot- the depth of the pedigree [2] because founders are as- spots that overlapped with signatures of selection [19]. sumed to be unrelated and non-inbred [3]. Conse- The majority of goat breeds analysed by Bertolini et al. quently, inbreeding produced by distant ancestors not [19] displayed low levels of inbreeding (FROH < 0.10), included in the pedigree is systematically ignored [4]. with the only exception of certain local breeds with Another disadvantage of quantifying inbreeding from small population sizes (e.g. Dutch Landrace goats) as pedigree data is that it provides bare expectations about well as of breeds with insular origins (e.g. Icelandic and the fraction of the genome which is identical-by-descent Malagasy goats) which happened to be highly inbred [19, (IBD) [3]. With the advent of high-density arrays of sin- 20]. Despite the fact that inbreeding depression can have gle nucleotide polymorphisms (SNPs), it has become adverse effects on the profitability of farmers and animal possible to estimate genomic inbreeding coefficients breeders, very few studies have investigated its conse- which circumvent these limitations [5]. Indeed, import- quences on goat production [22, 23]. In this regard, ant advantages of genomic inbreeding coefficients over Marete et al. [22] and Deroide et al. [23] estimated, with their genealogical counterparts are: (i) higher accuracy genealogical methods, the effect of inbreeding depres- to differentiate among individuals within the same pedi- sion on the production of Kenya Alpine and Murciano- gree, since variation due to Mendelian sampling is cap- Granadina goats, respectively, and they found that in tured [4], (ii) higher accuracy to quantify shared ancestry both populations such effect was negligible. of genetic haplotypes [4], and (iii) the ability to map in- The goals of the current work were: (i) to measure the breeding to specific genomic regions [6]. Different types levels of inbreeding in a Murciano-Granadina resource of genomic inbreeding coefficients have been imple- population by using different genomic coefficients, and mented. While inbreeding coefficients based on the (ii) to use this information to infer the impact of in- proportion of homozygous SNPs (FHOM) just reflect breeding depression on dairy phenotypes recorded in identity-by-state (IBS) allele-sharing proportions [6], co- this population. efficients (FROH) based on measuring the fraction of the genome covered by runs of homozygosity (ROH) esti- Methods mate IBD allele sharing [4–8], making possible to disen- Animal material and phenotyping tangle recent from ancient inbreeding [3–9]. The animal material comprised 1040 Murciano-Granadina The increase of inbreeding might have adverse con- female goats from 15 farms located in the autonomous re- sequences on the fitness of livestock populations due gion of Andalusia (Spain). Murciano-Ganadina is a local to the loss of genetic variability, which can entail a Spanish breed officially created in 1975 by the crossbreed- long-term reduction of genetic variance (due to the ing of Murciano and Granadina goats [24]. Currently, it has fixation of alleles) and, consequently, a slowing down a census of 115,105 heads (2020), and its remarkable adapt- of the rate of response to selection in breeding ability to harsh environments as well as its good milking schemes [10–12]. Moreover, incremented levels of in- performance (mean of 586 kg/lactation; 5.1% of fat and breeding might reduce the mean phenotypic perform- 3.6% of protein in milk) have made it a very popular breed ance of livestock populations, a phenomenon known in Spain and other countries (https://www.mapa.gob.es/es/ as inbreeding depression (reviewed by Leroy [2]). Al- ganaderia/ temas/zootecnia/razas-ganaderas/razas/cata- though inbreeding depression is particularly intense logo-razas/caprino/murciano-granadina/). for fitness and reproduction traits [11], there is broad Blood samples were extracted from goats using vac- evidence that it also decreases dairy and growth per- uum tubes coated with EDTA K3 anticoagulant and formances [2, 13–15]. Besides, susceptibility to certain stored at − 20 °C until processing. Phenotypic records diseases, such as mastitis, is increased in inbred ani- for milk yield and composition traits were recorded in mals [16, 17]. In Holstein cattle, a 1% increase of in- the framework of the selection program of the breeding is expected to cause a reduction of $22–24 Murciano-Granadina goat breed [24]. Only phenotypes of lifetime net income per individual [18], while in corresponding to the first parity (recorded between the sheep the average economic loss per ewe amounts to years 2009 and 2017) were taken into consideration. The $17 for moderate inbreeding and $36 when inbreed- following phenotypes were recorded in 817 goats: milk ing is close to 50% [18]. yield measured in kilograms at 210 days (MY210), 240 Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 3 of 14 days (MY240) and 305 days (MY305), the natural loga- only SNPs meeting the following requirements were rithm of the somatic cell count divided by 1000 (lnSCC, used in the downstream analyses: (i) mapping to auto- to convert this value into a somatic cell count please use somes, (ii) displaying a minor allele frequency of 0.05 or the formula: elnSCC × 103 cells/mL), fat percentage (FP), higher, (iii) not deviating significantly (P < 0.00001) from protein percentage (PP) and lactose percentage (LP). the Hardy-Weinberg equilibrium, and (iv) with a geno- Milk composition traits were standardised to a lactation type call rate over 98%. Data were filtered using PLINK of 210 days. Summary statistics of phenotypic records 1.9 [27]. In addition, a principal component analysis are displayed in Table 1. (PCA) was carried out with PLINK 1.9 [27] in order to assess population structure. The visualisation of such re- Generation of high throughput genotypic data sults was based on the first two components of the PCA. The isolation of genomic DNA was carried out following a salting-out protocol [25]. Three mL of whole blood Detection of runs of homozygosity were mixed with 4 volumes of Red Cell Lysis Solution The definition of ROH followed six criteria: (i) the mini- (Tris-HCl 10 mmol/L, pH = 6.5; EDTA 2 mmol/L; Tween mum length of ROH is 1 Mb, (ii) a ROH must contain 20 1%), and this mixture was centrifuged at 2000 × g. at least 15 SNPs, (iii) the density of SNPs per ROH was The supernatant was discarded and the pellet was resus- set to at least 1 SNP every 100 kb, (iv) the maximum dis- pended in 3 mL of lysis buffer (Tris-HCl 200 mmol/L, tance between consecutive SNPs is 250 Kb, (v) one het- pH = 8, EDTA 30 mmol/L, SDS 1%; NaCl 250 mmol/L) erozygous position per ROH is allowed, and (vi) one plus 100 μL proteinase K (20 mg/mL) and incubated for missing position per ROH is allowed. Said criteria were 3 h at 55 °C. The lysate was chilled, and 1 mL of ammo- established based on the density of the genotyping panel, nium acetate 10 mol/L was added to it. After centrifuga- with a mean distance between consecutive SNPs of tion at 2000 × g for 10 min, the supernatant (~ 4 mL) 51.73 Kb and a mean number of 19.35 SNPs/Mb (Add- was transferred to a new tube with 3 mL of isopropanol itional file 1: Table S1), and several of them are based on 96%, and this mixture was centrifuged at 2000 × g for 3 the ROH definition established by the AdaptMap Con- min. The resulting DNA pellet was washed with 3 mL of sortium [19, 20]. Runs of homozygosity were identified ethanol 70% and an additional centrifugation step at with the PLINK 1.9 software [27] using a sliding window 2000 × g for 1 min was performed. The DNA pellet was of 50 SNPs. dried at room temperature, and it was subsequently re- suspended in 1 mL of TE buffer (Tris-HCl 10 mmol/L, Analysis of the genomic patterns of homozygosity EDTA 1 mmol/L, pH = 8). Murciano-Granadina goats The proportion of homozygosity per site was estimated were genotyped with the Goat SNP50 BeadChip (Illu- as the ratio between the number of animals with homo- mina Inc., San Diego, CA, USA) by following the in- zygous genotypes for a particular SNP divided by the structions of the manufacturer. The goat ARS1 genome number of animals genotyped for that SNP. A sliding [26] was used as reference for inferring the genomic lo- window encompassing 30 SNPs was designed to esti- cation of the SNPs, and the position and the name of mate the average of this ratio, and chromosomal pat- each SNP were updated using the PLINK 1.9 software terns of homozygosity were visualised as Manhattan [27]. Only individuals with at least 95% of SNPs with plots using R [28]. genotype calls were taken into consideration. Moreover, Genomic coverage and distribution of ROH were also investigated. The patterns of ROH size and distribution Table 1 Summary statistics of seven milk production and along the genome were analysed and plotted using R composition traits recorded in 817 Murciano-Granadina goats [28]. Genomic regions in which ROH are prevalent, the Traitsa Mean SD so-called ROH hotspots, were identified by measuring MY210, kg 395.647 131.787 the proportion of animals that harbour a particular SNP MY240, kg 450.493 142.707 occurring within a ROH with regard to the total number MY305, kg 547.418 179.840 of animals genotyped for that SNP. Genomic regions containing the top 1% SNPs most commonly associated lnSCC 6.278 0.937 with ROH were classified as ROH hotspots [29, 30]. FP, % 5.190 0.766 Both highly homozygous regions and ROH hotspots PP, % 3.563 0.351 were compared, and overlapping segments were identi- LP, % 4.865 0.228 fied. Taking as reference the ARS1 goat genome [26], a MY210, milk yield at 210 days of lactation (kg); MY240, milk yield at 240 days genes mapping to these overlapping segments were ex- of lactation (kg); MY305, milk yield at 305 days of lactation (kg); lnSCC, natural tracted with the Biomart tool of Ensembl [31]. Database logarithm of the somatic cell count divided by 1000 (to convert this value into a somatic cell count use the formula: elnSCC × 103 cells/mL); FP, fat percentage; for Annotation, Visualization and Integration Discovery PP, protein percentage; LP, lactose percentage (DAVID) Bioinformatics Resources [32] and UniProt Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 4 of 14 [33] were used in Gene ontology analysis to identify estimated from the output of the --hardy command of over-represented (enriched) gene ontology (GO) terms PLINK 1.9 [27]. and KEGG pathways. Amongst other things, DAVID FYANG coefficient was estimated from the diagonal of provides biological context for long lists of genes by the matrix of genomic relationships of Yang based on assigning them to functionally related groups through the correlation between uniting gametes [37, 38]. It was the use of fuzzy clustering techniques [32]. To assess the calculated with the following formula: significance of gene-enrichment in annotation terms, a Fisher Exact is employed. Only terms with Max. Prob 1X s x2k −ð1 þ 2pk Þxk i þ 2p2k F YANGi ¼ ≤ 0.1, Min. Count ≥ 2 and P < 0.05 were considered as S k¼1 2pk ð1−pk Þ significantly enriched [32]. where, xk is the genotype of the individual i for the SNP Estimation of inbreeding coefficients k, and pk is the frequency of the reference allele in the Six inbreeding coefficients were estimated at the whole studied population. The command --ibc of PLINK 1.9 genome level based on genotypic data: [27] was used to estimate it. FHOM is defined as the proportion of SNPs with homo- The coefficients FHOM and FROH were also estimated zygous genotypes [6] and was estimated as: at the chromosomal level. H Oi F HOMi ¼ Inbreeding depression analyses S The effects of inbreeding depression on dairy traits were where HOi corresponds to the observed number of investigated using data from 817 goats with available homozygous genotypes for all the SNPs for each individ- phenotypes (MY210, MY240, MY305, lnSCC, FP, PP and ual i and S are the number of SNPs for which individual LP). Analyses were performed with the REMLF90 soft- i has genotype data. It was computed from the output of ware [39] to implement a restricted maximum likelihood the --het command of PLINK 1.9 software [27]. (REML) analysis approach, in which the phenotypic FROH was estimated as the proportion of the genome values of each trait in each individual are regressed onto covered by ROH by using the following formula: its inbreeding coefficient using a linear mixed model. P These analyses were performed to quantify inbreeding LROH i depression at the whole genome scale as well as on a F ROH i ¼ Lauto chromosome and regional basis. The model was fitted as follows: where LROHi corresponds to the sum of the lengths of all ROH present in each individual i, and Lauto is the total y ¼ Xβ þ Za þ e length of the autosomal goat genome covered by SNPs [5]. The same mathematical expression was used to calculate where y is the vector of observations for each phenotype, the genomic coverage of ROH with sizes smaller (FROHShort) β is a vector of fixed effects, including farm (15 levels), or larger (FROHLong) than 5 Mb. Such calculations were year of birth (10 levels), and the linear regression on F made to assess the relative importance of distant (FROHShort) as a covariate; a is the vector of additive genetic effects, versus recent (FROHLong) inbreeding [6], according to the e is the vector of random residual effects, and X and Z size of the genome and SNP density (Table 1). Assuming are incidence matrices relating fixed and random effects that the length of the ROH segments follows an exponen- to observations. tial distribution with a mean equal to 12 g , where g corre- The significance of the inbreeding effect on the ana- sponds to the number of generations to the closest lysed traits was determined by applying a two-tailed hy- common ancestor [34], and also assuming that goats have a pothesis test. A Z-statistic was estimated with the recombination rate of approximately 1 cM/Mb [35], FROH- following general formula: Short indicates the inbreeding of an individual from 10 to 50 generations in the past, while FROHLong estimates the in- −x −μ0 Z¼ breeding from 1 to 9 generations in the past. s:e: FIS coefficient of Wright [36] was calculated with the Here, x corresponds to the regression coefficient formula: representing the effect of inbreeding over each trait, and H Ei −H Oi μ0 is the coefficient of inbreeding corresponding to the F IS ¼ null hypothesis (in this case is equal to 0), and s.e. is the H Ei standard error. The transformation of Z-scores into P- Where H Ei and H Oi are the expected and observed values was accomplished with the function pnorm() im- heterozygosities of the individual i. This coefficient was plemented on R [28]. Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 5 of 14 The above analysis was performed by regressing each Sixty-six genes mapped to these regions (Additional file 6: phenotype onto four genomic inbreeding coefficients Table S4) and a functional enrichment analysis revealed (FHOM, FROH, FROHLong, and FROHShort). In order to detect 15 GO terms significantly enriched at the nominal level genomic regions associated with inbreeding depression, (P < 0.05). Particularly significant were GO terms related analyses at the chromosomal level were performed for with ferric and copper import into the cell (Add- traits that in the whole-genome analysis were identified itional file 7: Table S5). as significantly affected by inbreeding depression (P < 0.05). Following Saura et al. [6], inbreeding depres- Estimation of inbreeding coefficients sion was finely mapped by dividing chromosomes into Genomic inbreeding coefficients reached values of six segments and performing the analyses reported 0.601 ± 0.021 (FHOM); 0.053 ± 0.046 (FROH); 0.040 ± 0.041 above in each segment. As previously explained, genes (FROHLong); 0.014 ± 0.008 (FROHShort); − 0.016 ± 0.035 (FIS) mapping to genomic regions associated with inbreeding and 0.023 ± 0.047 (FYANG). The magnitude and disper- depression for a specific trait were retrieved using Bio- sion of these coefficients are shown in Fig. 4. The FHOM, mart [31], and their biological functions were assessed FROH, FROHLong, FIS and FYANG coefficients were highly with UniProt [33] and David Bioinformatics Resources correlated, being especially high the correlations between version 6.8 [32]. FHOM and FROH, FROH and FROHLong (r = 0.99, P< 2.2 × 10−16) and between FHOM and FIS (r = − 1, P < 2.2 × Results 10−16). In contrast, FROHShort showed the weakest corre- Assessment of homozygosity patterns in Murciano- lations with the remaining inbreeding coefficients (|r| = Granadina goats 0.33–0.64), although their statistical significance (P < A total of 46,689 SNPs and 1040 animals were selected 2.2 × 10−16) was very high (Table 2). to investigate the population structure and patterns of homozygosity of Murciano-Granadina goats. The PCA Measurement of inbreeding depression for dairy traits (Additional file 2: Fig. S1) indicated that goats clustered, The natural logarithmic transformation of the somatic to some extent, in accordance with their farm of origin. cell count divided by 1000 (lnSCC) as well as milk yield We detected 20,312 ROH that were classified as follows: at three different time points (MY210, MY240 and 11,325 had sizes of 0–5 Mb, while 5470 (5–10 Mb), 2695 MY305) showed significant evidence of inbreeding de- (10–20 Mb), 789 (20–50 Mb) and 33 (> 50 Mb) displayed pression when analysing FHOM and FROH. In contrast, no sizes above 5 Mb. The mean number of ROH per cat- significant effects of inbreeding were identified for coef- egory and per individual varied slightly across farms ficients based either on short or long ROH. Increases of (Fig. 1A and B). The mean ROH number was 19.53 ± 0.1 units of FHOM and FROH coefficients involved lnSCC 11.89 per individual, with an average length of 6.15 ± increments of 0.29 (P = 0.037) and 0.127 (P = 0.038) 2.05 Mb. As depicted in Fig. 2, the majority of the indi- units, respectively (Table 3). At the chromosomal level, viduals harboured less than 50 ROH, and ROH covered significant inbreeding depression for lnSCC was detected a small proportion of the genome (< 300 Mb, about 10% on CHI8 and CHI25 when regressed on either FHOM or of the genome). Only 2% of the individuals showed a FROH, while six additional chromosomes (CHI13, CHI14, genomic ROH coverage > 500 Mb, and 1% harboured CHI22, CHI24, CHI25 and CHI27) displayed inbreeding more than 50 ROH. Larger chromosomes encompassed depression for this trait exclusively when it was a greater number of ROH when compared to the smaller regressed onto FHOM (Additional file 8: Table S6). Four ones (Additional file 3: Fig. S2), and the correlation coef- regions containing 666 genes on chromosomes 8 (i.e. ficient between the number of ROH and chromosome CHI8:37,557,623–56,336,435 bp and CHI8:75,115,244– length was 0.92 (P < 0.05). 93,894,055 bp) and 25 (i.e. CHI25: 82,419–7,143,084 bp The genome-wide analysis of homozygosity, based on and CHI25:21,429,255–28,572,340 bp) displayed a sig- the proportion of homozygous individuals for each geno- nificant inbreeding depression for lnSCC based on FHOM typed position, made it possible to detect 38 genomic re- (Additional file 9: Table S7 and Additional file 10: Table gions with increased homozygosity (harbouring the top S8). After performing functional enrichment analysis, 1% of the most homozygous positions) that were scat- various GO terms and pathways involved in the immune tered on 20 goat chromosomes (CHI), i.e. CHI 1–8, 11, response were significantly enriched only at the nominal 13–18, 20, 21, 24, 26 and 29 (Fig. 3A and Add- level, including integrin domains (i.e. ITGAM, ITGAD, itional file 4: Table S2). Eight ROH hotspots mapping to ITGAX, ITGAL) and genes with protein kinase activity CHI 1, 2, 4, 6, 14, 16, and 17 were identified (Fig. 3B and chemotaxis and cell signalling activities (i.e. CCL19, Additional file 5: Table S3). One region (i.e. CHI4: CCL21, CCL24, CCL26, SYK and IL27). Several of these 42,552,375–48,378,207 bp) was consistently detected in genes display functions related with innate immunity the genome-wide analysis of homozygosity and ROH. and inflammatory response e.g. lymphocyte chemotaxis, Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 6 of 14 Fig. 1 A Number of ROH classified according to their length. Purple and yellow bars represent the counts of ROH shorter and longer than 5 Mb, respectively. B Number of ROH classified according to their length category by the farm of origin cellular response to interferon-γ, immunological synapse involved in mitochondrial and energetic processes, in- formation, positive regulation of chemotaxis, cellular re- cluding genes from the PPAR signalling pathway i.e. sponse to interleukin-1, monocyte chemotaxis etc. (Add- CYP27A1, ACADL, DBI and ACSL9 (Additional file 15: itional file 11: Table S9). Table S13). With regard to the milk yield traits, an increment of Although there is no substantial overlap between re- 0.1 units of FHOM or FROH involved a decrease of 48.50 gions associated with inbreeding depression and ROH kg (P = 0.011) and 20.49 kg (P = 0.014) of milk for hotspots, the CHI8:75–93.8 Mb region associated with MY210; 51.02 kg (P = 0.017) and 21.30 kg (P = 0.023) for inbreeding depression for lnSCC and the CHI2:45.5– MY240; and 55.72 kg (P = 0.032) and 23.20 kg (P = 0.040) 68.25 Mb region associated with inbreeding depression for MY305, respectively (Table 3). The analysis at the for milk yield show positional concordance with two of chromosomal level indicated significant inbreeding de- the 38 regions displaying high homozygosity (i.e. CHI8: pression for MY210 and MY240 on CHI2, CHI3 and 76.25–77.4 Mb and CHI2:56.13–57.17 Mb). CHI11, while for MY305 inbreeding depression was sig- nificant only on CHI2 (Additional file 12: Table S10). Discussion The region on CHI2, comprised between 22,751,824 and Low inbreeding in Murciano-Granadina goats 68,255,473 bp displayed significant inbreeding for all For the majority of Murciano-Granadina goats, ROH three milk yield traits (Additional file 13: Table S11). A number and total length were below 50 and 350 Mb, re- total of 355 genes mapped to this genomic region (Add- spectively. Moreover, short ROH (< 5 Mb) were more itional file 14: Table S12). The functional enrichment abundant than the medium or long ones (Fig. 1A). These analysis evidenced an overrepresentation of genes patterns are pretty consistent with what has been Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 7 of 14 Fig. 2 Relationship between ROH number and length in 1040 Murciano-Granadina goats genotyped with the Goat SNP50 BeadChip observed by Bertolini et al. [19] in a worldwide sample breed. Widespread use of artificial insemination in re- of goat breeds. Indeed, Bertolini et al. [19] showed that productive management and intensive selection were goat breeds from Southern Europe had, on average, 49 implemented in the Murciano-Granadina breed a few ROH per individual while the genomic coverage per in- decades ago [24], so their impact on genetic diversity dividual was 183.47 Mb. In contrast, breeds from North- and inbreeding has probably been quite limited so far. ern Europe showed higher levels of homozygosity with While correlations between FHOM, FROH, FROHLong, FIS 98 ROH per individual and genomic coverage of 479.17 and FYANG were high, FROHShort displayed the lowest cor- Mb. relations with the remaining inbreeding coefficients, in The inbreeding coefficients FROH, FROHShort, and FROH- line with previous studies focused on other livestock Long of Murciano-Granadina goats were mainly in the species [6, 13, 41]. Short ROH are mainly derived from range of 0 to 0.05. In their study, Bertolini et al. [19] re- ancient inbreeding events [5, 42] and do not reflect the ported that about 60% of a worldwide sample of goat whole autozygosity of the sample. It is also possible that breeds displayed low FROH coefficients (< 0.10), while the several of these homozygous tracks are identical by state remaining ~ 30% and ~ 10% of breeds showed moderate and not by descent, being produced by a low recombin- (0.10 < FROH < 0.20) or high (> 0.20) FROH values. The ation rate or high linkage disequilibrium in unrelated an- patterns of low homozygosity that we have observed in cestors [30]. Besides, when working with medium Murciano-Granadina goats contrast strongly with what density genotype arrays (e.g. 50 K SNPs) the detection of has been reported in certain local breeds, such as short ROH can become quite difficult [43], thus limiting Mallorquina, Pyrenean, and Valdostana, which have the ability to infer the true proportion of short vs. long undergone sharp population bottlenecks [19, 20]. Low ROH in the genome. inbreeding and homozygosity in the Murciano- Granadina breed are probably explained by its large Several ROH hotspots are detected in the genomes of census size (> 100,000 individuals registered in the herd- Murciano-Granadina goats book), the absence of population bottlenecks, and its We have identified several genomic regions in which broad geographic distribution encompassing more than ROH are particularly frequent (ROH hotspots). Similar 4000 farms across Spain (https://www.mapa.gob.es/es/ patterns were found when the proportion of homozy- ganaderia/temas/zootecnia/razas-ganaderas/razas) and gous individuals was analysed at a genome-wide level i.e. other countries. Noteworthy, the Murciano-Granadina 35.46% of the most common homozygous regions over- breed was founded by crossing, during the 1970s, two lapped with ROH hotspots and these overlapping re- Murciana and Granadina populations with different gions represented 35.42% of the total ROH hotspots. No historical origins [40]. Although genetic differentiation positional coincidence was detected between ROH hot- between these two populations was weak [40], this spots identified by us and those reported by Bertolini admixture event probably contributed to increase the et al. [19] in a worldwide sample of goat populations. heterozygosity of the resulting Murciano-Granadina This finding agrees well with what has been previously Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 8 of 14 Fig. 3 A Proportion of individuals with homozygous genotypes for each SNP marker. The y-axis displays the proportion of individuals for which a specific SNP displays a homozygous genotype, while the x-axis corresponds to the positional coordinates of SNPs distributed in the 29 caprine autosomes. B ROH hotspots identified in the population of Murciano-Granadina goats under study. The y-axis displays the frequency at which a given SNP is found within a ROH in the population; while the x-axis corresponds to the positional coordinates of SNPs distributed in the 29 caprine autosomes. Markers above the red line are in the top 1% of each category (homozygosity or frequency of being within a ROH). Markers highlighted in green are located in genomic regions consistently identified as regions of high homozygosity and ROH hotspots observed in sheep [44]. Indeed, ROH hotspots are ROH hotspots contain genes with diverse functions produced by factors such as positive selection and inver- Regarding the gene content of genomic segments co- sions suppressing recombination, that can differ sub- localising with both high homozygosity regions and stantially from population to population [42]. For ROH hotspots (CHI4:42,552,375–48,378,207 bp), the instance, the characterisation of the products of 5860 functional enrichment analysis highlighted several gene female meioses in Drosophila melanogaster by geno- ontology terms with nominally significant enrichment typing more than 100 million SNPs made it possible (P-value < 0.05) (Additional file 7: Table S5). From this to detect 106,964 recombination events displaying a list of genes, it is worth emphasising STEAP1, STEAP2 remarkable degree of intra-specific variation [45]. Fac- and STEAP4 that encode homonymous metalloreduc- tors such as GC content, gene density, distribution of tases which facilitate the cellular uptake of iron and simple repeats and transposable elements, structural copper [47]. These proteins modulate the effects of variation, and the presence of diverse poorly- intracellular oxidative stress and inflammation and are characterised sequence motifs might explain the re- involved in multiple biological pathways related with gional variation of the recombination rate across indi- molecular trafficking in the endocytic and exocytic path- viduals and populations [46]. ways, metabolism, control of cell proliferation and Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 9 of 14 Fig. 4 Boxplots depicting the magnitude and dispersion of molecular inbreeding FROH, FROHLong, FROHShort, FHOM, FIS and FYANG coefficients estimated in 1040 female Murciano-Granadina goats. Differences in magnitude between FHOM and the other molecular coefficients are due to the fact that they indicate identity-by-state and identity-by-descent allele-sharing proportions, respectively apoptosis and tumour progression [48]. We also found apocrine nature of milk secretion in goats increases the genes with metabolic functions such as insulin-like proportion of cytoplasmatic particles in milk, a feature, growth factor binding protein 3 (IGFBP3) and insulin- that depending on the measurement method of choice, like growth factor binding protein 1 (IGFBP1) that par- might increase the somatic cell count [52]. Besides, som- ticipate in the growth and postnatal development of cat- atic cell count is modulated by many factors including tle [49, 50]. These results are quite concordant with the occurrence of bacterial infections, stress, oestrous those reported by Mastrangelo and colleagues [30], who cycle phase, diet etc. [53]. According to our results showed that ROH islands identified in Italian bovine (Table 3), inbreeding depression increased lnSCC, a fea- breeds contained genes with heterogeneous functions re- ture that is considered adverse because high lnSCC lated to milk production, reproduction, immune re- values are often associated with subclinical and clinical sponse, and resistance/susceptibility to infection and mastitis [51, 52]. diseases. The magnitude of inbreeding depression for lnSCC estimated from FROHShort strongly differed from estimates Effect of inbreeding depression on milk traits based on the other inbreeding coefficients (Table 3), a The lnSCC in the Murciano-Granadina population finding consistent with the moderate (|r| = 0.33–0.64) cor- under study averaged 6.25 ± 0.93 units, which is higher relations between FROHShort and other molecular inbreed- than values reported in primiparous goats from the Al- ing coefficients (Table 2). Noteworthy, recent rather than pine (5.09 ± 1.36 units) and Saanen (5.32 ± 1.19 units) ancestral inbreeding is the main cause of inbreeding de- breeds [51]. Compared with cows and ewes, goats dis- pression in mammalian populations [54, 55]. Besides, long play higher numbers of somatic cells in milk. Indeed, the stretches of homozygosity usually contain a higher proportion of deleterious mutations than the shorter ones Table 2 Pearson correlations between molecular inbreeding because they are more recent, so deleterious variation has coefficients (F) estimated in 1040 Murciano-Granadina goats not been yet purged by purifying selection [56]. F coefficient FHOM FROH FROHLong FROHShort FIS FYANG In Murciano-Granadina goats, a previous study per- FHOM 1 formed by Deroide et al. [23] reported a low percentage of inbreeding (average F = 0.24%). Milk production FROH 0.99 1 showed a positive quadratic correlation with inbreeding FROHLong 0.97 0.99 1 levels, but such effect was not significant. Deroide et al. FROHShort 0.64 0.59 0.45 1 [23] also reported that milk fat and dry extract contents FIS –1 – 0.99 – 0.97 – 0.64 1 experienced a slight increase due to inbreeding. In our FYANG 0.86 0.89 0.92 0.33 – 0.86 1 study, the dairy trait mostly affected by inbreeding de- *All correlation coefficients (lower part of the matrix) were significant pression was lnSCC. Consistently, Doekes et al. [13] re- (P < 2.2 × 10−6) ported that a 1% increase of FROH involved a 0.86 ± 0.28 Table 3 Inbreeding depression estimates for milk yield and composition traits expressed as the change of the phenotypic mean per 0.1 units increase of the corresponding inbreeding coefficient (95% confidence intervals are displayed in brackets) Traitb FHOM FROH FROHLong FROHShort Mean ± s.e.a P-value Mean ± s.e.a P-value Mean ± s.e.a P-value Mean ± s.e.a P-value (CI 95%) (CI 95%) (CI 95%) (CI 95%) Luigi-Sierra et al. Journal of Animal Science and Biotechnology MY210 – 48.503 ± 21.124 (− 89.906 to – 7.101) 0.011 – 20.492 to ±9.356 (− 38.830 to – 2.154) 0.014 21.231 ± 25.901 (− 29.536 to 71.998) 0.206 19.751 ± 83.308 (− 143.533 to 183.034) 0.406 MY240 – 51.017 ± 24.159 (− 98.369 to – 3.666) 0.017 – 21.296 to ±10.702 (− 42.272 to – 0.320) 0.023 34.358 ± 30.841 (− 26.090 to 94.806) 0.133 10.158 ± 93.948 (− 173.980 to 194.296) 0.457 MY305 – 55.719 ± 29.987 (− 114.493 to 3.055) 0.032 – 23.196 to ±13.284 (− 49.233 to 2.841) 0.040 48.052 ± 47.583 (− 45.211 to 141.315) 0.156 – 16.834 ± 118.042 (− 248.196 to 214.527) 0.443 lnSCC 0.290 ± 0.162 (− 0.027 to 0.608) 0.037 0.127 ± 0.072 (− 0.013 to 0.268) 0.038 – 0.268 ± 0.635 (− 1.514 to 0.977) 0.336 – 0.368 ± 0.649 (− 1.640 to 0.905) 0.286 (2022) 13:35 FP 0.119 ± 0.121 (− 0.117 to 0.356) 0.161 0.046 ± 0.053 (− 0.059 to 0.150) 0.196 – 0.016 ± 0.209 (− 0.426 to 0.394) 0.469 – 0.031 ± 0.428 (− 0.871 to 0.809) 0.471 LP 0.029 ± 0.044 (− 0.057 to 0.116) 0.252 0.015 ± 0.019 (− 0.024 to 0.053) 0.226 – 0.090 ± 0.172 (− 0.428 to 0.247) 0.300 – 0.073 ± 0.169 (− 0.404 to 0.257) 0.332 PP 0.071 ± 0.068 (− 0.062 to 0.204) 0.147 0.035 ± 0.030 (− 0.024 to 0.093) 0.122 0.039 ± 0.106 (− 0.169 to 0.246) 0.358 – 0.157 ± 0.257 (− 0.660 to 0.346) 0.270 a s.e. standard error; CI 95, 95% confidence interval b MY210; milk yield at 210 days of lactation (kg); MY240; milk yield at 240 days of lactation (kg); MY305; milk yield at 305 days of lactation (kg); lnSCC, natural logarithm of the milk somatic cell count divided by 1000 (to convert this value into a somatic cell count use the formula: elnSCC × 103 cells/mL); FP, milk fat percentage (%); LP, milk lactose percentage (%); PP, milk protein percentage (%) * numbers in bold correspond to P-values below 0.05 Page 10 of 14 Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 11 of 14 unit increase in somatic cell score (days 150 through to essential for the migration and adhesion of polymorpho- 400) recorded in Dutch Holstein-Friesian dairy cattle. In nuclear leukocytes to infection sites [68]. Iranian cattle, individuals with high inbreeding coeffi- As previously indicated, genes related to chemotaxis cients tended to have higher somatic cell scores than an- were significantly enriched at the nominal level (P-value imals with low inbreeding coefficients [57], and similar < 0.05) (Additional file 11: Table S9). This functional results have been reported for Canadian Holstein cattle category is mainly represented by chemokine genes, such [58]. Doekes et al. [13] indicated that ancient inbreeding as chemokine ligand 19 (CCL19), 21 (CCL21), 24 was the main contributor to inbreeding depression for (CCL24), 26 (CCL26) and 27 (CCL27). Chemokines are somatic cell score, although such effect was not signifi- essential for the development of the innate immune re- cant. Somatic cell score is an indicator of the health sta- sponse since they orchestrate and control the migration tus of the mammary gland and substantial increases are of the immune cells (macrophages, monocytes, neutro- observed in individuals suffering from mastitis [59]. In- phils, etc.) to sites of infection [69]. During the first breeding has been reported to significantly reduce resist- stages of mastitis, chemokines contribute to the stimula- ance against pathogens in multiple organisms [60–62], tion of the cellular immune response against the invad- so the significant inbreeding depression observed for ing pathogen until acute-phase proteins are expressed lnSCC in Murciano-Granadina goats might be explained, [70]. Marsland et al. [71] described how chemokines at least in part, by the weakening of the immune de- CCL19 and CCL21 participate in the maturation of den- fences of the mammary gland. Thus, homozygosity for dritic cells, allowing them to leverage the T cell re- deleterious mutations might result in the partial or total sponse. These chemokines also enhance the migration of inactivation of genes related with immunity, and low leukocytes through lymph and blood circulation and variability might also compromise the effectiveness of stimulate the production of pathogen-induced proin- the immune response [63]. flammatory cytokines [71]. In the gene set associated with inflammatory response (P-value = 0.03), we detected the interleukin 27 (IL27) gene which encodes a molecule Genomic regions associated with inbreeding depression with both pro and anti-inflammatory effects, thus enhan- for lnSCC contain several genes related with immunity cing the immune response and, at the same time, pre- When we investigated which enriched clusters are present venting tissue damage caused by inflammation [72]. in the set of genes mapping to chromosomes 8 (37–56 Mb Infections caused by Gram-negative bacteria induce IL27 and 75–93 Mb) and 25 (0.082–7 Mb and 21–28 Mb) re- production, and this cytokine interacts with monocytes gions associated with inbreeding depression for lnSCC increasing TLR4 expression and enhancing the LPS- (Additional file 10: Table S8), we found several genes induced inflammatory response [73]. Moreover, IL27 has assigned to gene ontologies highly connected with immun- an autocrine effect on macrophages and monocytes ity, e.g. integrin-mediated signalling pathway, lymphocyte resulting in the amplification of the inflammatory chemotaxis, monocyte chemotaxis, immunological synapse response via cytokine secretion [72]. formation, chemokine activity, etc. (Additional file 11: Table The genomic region displaying inbreeding depres- S9). The spleen tyrosine kinase (SYK) protein forms part of sion for milk yield (MY210, MY240 and MY305) was the integrins cluster and maps to CHI8: 86,755,291– significantly enriched at the nominal level (P-value < 86,861,895 bp (Additional file 11: Table S9). One of the 0.05) with genes associated with multiple unrelated functions of the SYK molecule is to stimulate the biological processes (Additional file 15: Table S13). phosphorylation of Toll-like receptor 4 (TLR4) [64], which The PPAR signalling pathway (P-value = 0.032) influ- recognises bacterial lipopolysaccharide and induces inflam- ences milk production and composition in cattle. Bai matory and immune responses [65]. This gene has been and collaborators [74] described an overrepresentation described as highly variable in cattle [66], and many TLR4 of genes from the PPAR signalling pathway in cows polymorphisms and haplotypes have been associated with with high daily milk yield in comparison with low milk somatic cell count and susceptibility to mastitis [67]. yielders. Besides, genes from this pathway are upregu- Moreover, the SYK protein induces the recognition of path- lated in cows in the lactation peak when compared to ogens and cell adhesion and platelet activation [65], and it those in the dry period [75], suggesting a role of also affects the proliferation of mammary epithelial cells at these genes not only in the determinism of fat com- several stages of the milking cycle [66]. In the same position [76] but also of milk yield. enriched gene ontologies, we have detected the integrin subunit αM gene (ITGAM also known as CD11b) which Conclusions maps to CHI25: 27,221,164–27,264,350 bp and encodes a Murciano-Granadina goats display low levels of in- receptor for lipopolysaccharide [68]. Signalling mediated by breeding (mean FROH = 0.053 ± 0.046), a finding con- TLR4 activates the synthesis of ITGAM/CD11b, which is sistent with the large census size and demographic Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 12 of 14 history of this breed. Four genomic regions on CHI8: regions associated with inbreeding depression for the natural logarithm 37,557,623–56,336,435, CHI8:75,115,244–93,894,055 of the somatic cell count divided by 1000 (lnSCC). bp, CHI25:82,419–7,143,084 and CHI25:21,429,255– Additional file 11: Table S9. Functional enrichment analysis of genes 28,572,340 bp were associated with inbreeding depres- from the goat chromosome 8 (37–56 Mb and 75–93 Mb) and 25 (0.082–7 Mb and 21–28 Mb) regions associated with inbreeding depression for the sion for lnSCC. Moreover, one region on CHI2: natural logarithm of the somatic cell count divided by 1000 (lnSCC). 22,751,824–68,255,473 was consistently associated Additional file 12: Table S10. Inbreeding depression estimates (s.e: with inbreeding depression for three milk yield traits standard error) per goat chromosome (CHI) for milk yield at 210, 240 and (MY210, MY240 and MY310). Genes encoding integ- 305 days (MY210, MY240 and MY305, measured in kg) and their 95% confidence intervals (C.I). rins, chemokines and pathogen recognition receptors, Additional file 13: Table S11. Inbreeding depression estimates (s.e: which play relevant roles in the elicitation of innate standard error) for milk yield at 210, 240 and 305 days (MY210, MY240 immune responses against microbes, mapped to and MY305, measured in kg) in specific regions of goat chromosomes regions associated with inbreeding depression for (CHI) 2 and their 95% confidence intervals (C.I). lnSCC. These results suggest that keeping inbreeding Additional file 14: Table S12. Genes mapping to goat chromosome (CHI) 2: 22–68 Mb region associated with inbreeding depression for milk to a minimum, through an adequate reproductive yield at 210, 240 and 305 days (MY210, MY240 and MY305). management, might be a useful approach to decrease Additional file 15: Table S13. Functional enrichment analysis of genes the incidence of mastitis in Murciano-Granadina from the goat chromosome (CHI) 2: 22–68 Mb region associated with goats. inbreeding depression for milk yield at 210, 240 and 305 days (MY210, MY240 and MY305). Abbreviations Acknowledgements IBD: Identical-by-descent; SNPs: Single nucleotide polymorphisms; Many thanks to CAPRIGRAN for carrying out phenotype recording and blood IBS: Identity-by-state; ROH: Runs of homozygosity; CHI: Chromosomes; FP: Fat sample collection in Murciano-Granadina goats. percentage; PP: Protein percentage; LP: Lactose percentage; PCA: Principal component analysis; DAVID: Annotation Visualization and Integration Authors’ contributions Discovery; GO: Gene ontology; REML: Restricted maximum likelihood; MA, JJ, JVD and VL designed the study. JFA, JVD and AM coordinated all SYK: Spleen tyrosine kinase tasks involved in phenotype recording. VL did all DNA extractions. MGLS, MS and AF carried out homozygosity and inbreeding depression analyses, with the cooperation of DG. MGLS and MA wrote the first draft of the paper with Supplementary Information the cooperation of MS and AF. All authors read and approved the content of The online version contains supplementary material available at https://doi. the paper. org/10.1186/s40104-022-00684-5. Funding This research was funded by the European Regional Development Fund Additional file 1: Table S1. Number and density of SNPs per (FEDER)/Ministerio de Ciencia e Innovación - Agencia Estatal de chromosome in a population of 1040 Murciano-Granadina goats geno- Investigación/Project Reference grant: PID2019-105805RB-I00 and by the typed with the Goat SNP50 BeadChip. CERCA Programme/Generalitat de Catalunya. We also acknowledge the sup- Additional file 2: Figure S1. Principal component analysis of 1040 port of the Spanish Ministerio de Ciencia e Innovación for the Center of Ex- Murciano-Granadina goats distributed in 15 farms (each farm is indicated cellence Severo Ochoa 2020–2023 (CEX2019–000902-S) grant awarded to the with a different colour). Centre for Research in Agricultural Genomics (CRAG, Bellaterra, Spain). We Additional file 3: Figure S2. Number of ROH per chromosome also acknowledge the support of the CERCA programme of the Generalitat (represented as yellow bars, left axis) and the percentage of each de Catalunya. Dailu Guan was funded by a PhD fellowship from the China chromosome covered by ROH (represented by a red line, right axis) in Scholarship Council (CSC). Maria Luigi-Sierra was funded with a PhD fellow- 1040 Murciano-Granadina goats. ship Formación de Personal Investigador (BES-C-2017-079709) awarded by the Spanish Ministry of Economy and Competitivity. Additional file 4: Table S2. Genomic regions associated with the top 1% of homozygosity for each SNP marker in a population of 1040 Availability of data and materials Murciano-Granadina goats genotyped with the Goat SNP50 BeadChip. Goat SNP50 BeadChip genotypes and milk production phenotypes are Additional file 5: Table S3. ROH hotspots identified in the genomes of accessible at Figshare (https://doi.org/10.6084/m9.figshare.18095825). 1040 Murciano-Granadina goats. Additional file 6: Table S4. Genes mapping to regions consistently Declarations identified as ROH hotspots and regions with high homozygosity. Ethics approval and consent to participate Additional file 7: Table S5. Functional enrichment analysis of genes Blood collection is a routine procedure performed by CAPRIGRAN, so it does mapping to regions consistently identified as ROH hotspots and regions not require approval by the Ethics Committee on Animal and Human with high homozygosity. Experimentation of the Universitat Autònoma de Barcelona. Additional file 8: Table S6. Inbreeding depression estimates (s.e: standard error) per goat chromosome (CHI) for the natural logarithm of Consent for publication the somatic cell count divided by 1000 (lnSCC) and their 95% confidence Not applicable. intervals (C.I). Additional file 9: Table S7. Inbreeding depression estimates (s.e: Competing interests standard error) for the natural logarithm of the somatic cell count The authors declare that they have no competing interests. divided by 1000 (lnSCC) in specific regions of goat chromosomes (CHI) 8 and 25 and their 95% confidence intervals (C.I). Author details 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Additional file 10: Table S8. Genes mapping to goat chromosome (CHI) 8 (37–56 Mb and 75–93 Mb) and 25 (0.082–7 Mb and 21–28 Mb) Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain. 2Departamento de Mejora Genética Animal, INIA, Carretera de la Coruña km 7,5, 28040 Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 13 of 14 Madrid, Spain. 3Departamento de Genética, Universidad de Córdoba, 14071 population history and adaptation of goats. Genet Sel Evol. 2018;50(1):59. Córdoba, Spain. 4Department of Veterinary Medicine, University of Bari ‘‘Aldo https://doi.org/10.1186/s12711-018-0424-8. Moro”, 62 per Casamassima km. 3, 70010 Valenzano, SP, Italy. 5Group of 20. Cardoso TF, Amills M, Bertolini F, Rothschild M, Marras G, Boink G, et al. Research in Ruminants (G2R), Department of Animal and Food Science, Patterns of homozygosity in insular and continental goat breeds. Genet Sel Universitat Autònoma de Barcelona (UAB), Bellaterra, Barcelona, Spain. Evol. 2018;50(1):56. https://doi.org/10.1186/s12711-018-0425-7. 6 Departament de Ciència Animal i dels Aliments, Facultat de Veterinària, 21. Nandolo W, Mészáros G, Banda LJ, Gondwe TN, Lamuno D, Mulindwa HA, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain. et al. Timing and extent of inbreeding in African goats. Front Genet. 2019; 10:537. https://doi.org/10.3389/fgene.2019.00537. Received: 3 August 2021 Accepted: 19 January 2022 22. Marete AG, Jung’a JO, Mosi RO. An analysis of inbreeding levels and factors affecting growth and reproductive traits in the Kenya Alpine dairy goat. Livest Res Rural Dev. 2011;23:259. 23. Deroide CAS, Jacopini LA, Delgado JV, Léon JM, Brasil LHA, Ribeiro MN. References Inbreeding depression and environmental effect on milk traits of the 1. Lush JL. The consequences and measurement of inbreeding. In: Animal Murciano-Granadina goat breed. Small Rumin Res. 2016;134:44–8. https:// breeding plans. Ames: The Iowa State College Press; 1943. 259–288. doi.org/10.1016/j.smallrumres.2015.12.008. 2. Leroy G. Inbreeding depression in livestock species: review and meta- 24. Delgado JV, Landi V, Barba CJ, Fernández J, Gómez MM, Camacho ME, et al. analysis. Anim Genet. 2014;45(5):618–28. https://doi.org/10.1111/age.12178. Murciano-Granadina goat: a Spanish local breed ready for the challenges of 3. Kardos M, Luikart G, Allendorf FW. Measuring individual inbreeding in the the twenty-first century. In: Simões J, Gutiérrez C, editors. Sustainable goat age of genomics: marker-based measures are better than pedigrees. production in adverse environments: volume II. Springer: Cham; 2018. p. Heredity. 2015;115(1):63–72. https://doi.org/10.1038/hdy.2015.17. 205–19. https://doi.org/10.1007/978-3-319-71294-9_15. 4. Keller MC, Visscher PM, Goddard ME. Quantification of inbreeding due to 25. Miller SA, Dykes DD, Polesky HF. A simple salting out procedure for distant ancestors and its detection using dense single nucleotide extracting DNA from human nucleated cells. Nucleic Acids Res. 1988;16(3): polymorphism data. Genetics. 2011;189(1):237–49. https://doi.org/10.1534/ 1215. https://doi.org/10.1093/nar/16.3.1215. genetics.111.130922. 26. Bickhart DM, Rosen BD, Koren S, Sayre BL, Hastie AR, Chan S, et al. Single- 5. McQuillan R, Leutenegger AL, Abdel-Rahman R, Franklin CS, Pericic M, molecule sequencing and chromatin conformation capture enable de novo Barac-Lauc L, et al. Runs of homozygosity in European populations. Am J reference assembly of the domestic goat genome. Nat Genet. 2017;49(4): Hum Genet. 2008;83(3):359–72. https://doi.org/10.1016/j.ajhg.2008.08.007. 643–50. https://doi.org/10.1038/ng.3802. 6. Saura M, Fernández AIA, Varona L, Fernández AIA, De Cara MÁR, Barragán C, et al. Detecting inbreeding depression for reproductive traits in Iberian pigs 27. Chang CC, Chow CC, Tellier LCAM, Vattikuti S, Purcell SM, Lee JJ. Second- using genome-wide data. Genet Sel Evol. 2015;47(1):1. https://doi.org/10.11 generation PLINK: rising to the challenge of larger and richer datasets. 86/s12711-014-0081-5. Gigascience. 2015;4(1):7. https://doi.org/10.1186/s13742-015-0047-8. 7. Broman KW, Weber JL. Long homozygous chromosomal segments in 28. R Core Team. R: A language and environment for statistical computing. reference families from the Centre d’Etude du Polymorphisme Humain. Am Vienna, Austria: R Foundation for Statistical Computing. ISBN 3–900051–07- J Hum Genet. 1999;65(6):1493–500. https://doi.org/10.1086/302661. 0, URL http://www.R-project.org/ 8. Pemberton TJ, Absher D, Feldman MW, Myers RM, Rosenberg NA, Li JZ. 29. Purfield DC, McParland S, Wall E, Berry DP. The distribution of runs of Genomic patterns of homozygosity in worldwide human populations. Am J homozygosity and selection signatures in six commercial meat sheep Hum Genet. 2012;91(2):275–92. https://doi.org/10.1016/j.ajhg.2012.06.014. breeds. PLoS ONE. 2017;12(5):e0176780. https://doi.org/10.1371/journal. 9. Purfield DC, Berry DP, McParland S, Bradley DG. Runs of homozygosity and pone.0176780. population history in cattle. BMC Genet. 2012;13(1):70. https://doi.org/10.11 30. Mastrangelo S, Tolone M, Sardina MT, Sottile G, Sutera AM, Di Gerlando R, 86/1471-2156-13-70. et al. Genome-wide scan for runs of homozygosity identifies potential 10. Kristensen TN, Sørensen AC. Inbreeding - lessons from animal breeding, candidate genes associated with local adaptation in Valle del Belice sheep. evolutionary biology and conservation genetics. Anim Sci. 2005;80(2):121– Genet Sel Evol. 2017;49(1):84. https://doi.org/10.1186/s12711-017-0360-z. 33. https://doi.org/10.1079/ASC41960121. 31. Kinsella RJ, Kähäri A, Haider S, Zamora J, Proctor G, Spudich G, et al. 11. Charlesworth D, Willis JH. The genetics of inbreeding depression. Nat Rev Ensembl BioMarts: A hub for data retrieval across taxonomic space. Genet. 2009;10(11):783–96. https://doi.org/10.1038/nrg2664. Database. 2011;2011:bar030. 12. Howard JT, Pryce JE, Baes C, Maltecca C. Invited review: inbreeding in the 32. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of genomics era: inbreeding, inbreeding depression, and management of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1): genomic variability. J Dairy Sci. 2017;100(8):6009–24. https://doi.org/10.3168/ 44–57. https://doi.org/10.1038/nprot.2008.211. jds.2017-12787. 33. Bateman A. UniProt: a worldwide hub of protein knowledge. Nucleic Acids 13. Doekes HP, Veerkamp RF, Bijma P, De Jong G, Hiemstra SJ, Windig JJ. Res. 2019;47(D1):D506–15. https://doi.org/10.1093/nar/gky1049. Inbreeding depression due to recent and ancient inbreeding in Dutch 34. Thompson EA. Identity by descent: variation in meiosis, across genomes, Holstein-Friesian dairy cattle. Genet Sel Evol. 2019;52(1):64. https://doi.org/1 and in populations. Genetics. 2013;194(2):301–26. https://doi.org/10.1534/ 0.1186/s12711-020-00583-1. genetics.112.148825. 14. Pryce JE, Haile-Mariam M, Goddard ME, Hayes BJ. Identification of genomic 35. Dumont BL, Payseur BA. Evolution of the genomic rate of recombination in regions associated with inbreeding depression in Holstein and Jersey dairy mammals. Evolution. 2008;62(2):276–94. https://doi.org/10.1111/j.1558-564 cattle. Genet Sel Evol. 2014;46(1):71. https://doi.org/10.1186/s12711-014- 6.2007.00278.x. 0071-7. 36. Wright S. The genetical structure of populations. Ann Eugenics. 1951;15(4): 15. Thompson JR, Everett RW, Hammerschmidt NL. Effects of inbreeding on 323–54. https://doi.org/10.1111/j.1469-1809.1949.tb02451.x. production and survival in Holsteins. J Dairy Sci. 2000;83(9):2131–8. https:// 37. Villanueva B, Fernández A, Saura M, Caballero A, Fernández J, Morales- doi.org/10.3168/jds.S0022-0302(00)75096-X. González E, et al. The value of genomic relationship matrices to estimate 16. Sørensen AC, Madsen P, Sørensen MK, Berg P. Udder health shows levels of inbreeding. Genet Sel Evol. 2021;53(1):42. https://doi.org/10.1186/ inbreeding depression in Danish Holsteins. J Dairy Sci. 2006;89(10):4077–82. s12711-021-00635-0. https://doi.org/10.3168/jds.S0022-0302(06)72452-3. 38. Yang J, Benyamin B, McEvoy BP, Gordon S, Henders AK, Nyholt DR, et al. 17. Croquet C, Mayeres P, Gillon A, Vanderick S, Gengler N. Inbreeding Common SNPs explain a large proportion of the heritability for human depression for global and partial economic indexes, production, type, and height. Nat Genet. 2010;42(7):565–9. https://doi.org/10.1038/ng.608. functional traits. J Dairy Sci. 2006;89(6):2257–67. https://doi.org/10.3168/jds. 39. Misztal I, Tsuruta S, Strabel T, Auvray B, Druet T, Lee DH. BLUPF90 and S0022-0302(06)72297-4. related programs (BGF90). In Proceedings of the 7th World Congress on 18. Ercanbrack SK, Knight AD. Effects of inbreeding on reproduction and wool Genetics Applied to Livestock Production: 19-23 August 2002; Montpellier. production of Rambouillet, Targhee, and Columbia ewes. J Anim Sci. 1991; 2002:28.07. 69(12):4734–44. https://doi.org/10.2527/1991.69124734x. 40. Martínez A, Rocha L, Quiroz J, Delgado BJ. Estudio de la diversidad genética 19. Bertolini F, Cardoso TF, Marras G, Nicolazzi EL, Rothschild MF, Amills M. intrarracial de la cabra Murciano-Granadina con microsatélites de ADN. Arch Genome-wide patterns of homozygosity provide clues about the Zootec. 2007;56:417–20. Luigi-Sierra et al. Journal of Animal Science and Biotechnology (2022) 13:35 Page 14 of 14 41. Biscarini F, Mastrangelo S, Catillo G, Senczuk G, Ciampolini R. Insights into 62. Smallbone W, van Oosterhout C, Cable J. The effects of inbreeding on genetic diversity, runs of homozygosity and heterozygosity-rich regions in disease susceptibility: Gyrodactylus turnbulli infection of guppies, Poecilia Maremmana semi-feral cattle using pedigree and genomic data. Animals. reticulata. Exp Parasitol. 2016;167:32–7. https://doi.org/10.1016/j.exppara.201 2020;10(12):2285. https://doi.org/10.3390/ani10122285. 6.04.018. 42. Ceballos FC, Joshi PK, Clark DW, Ramsay M, Wilson JF. Runs of 63. Siljestam M, Rueffler C. Heterozygote advantage can explain the homozygosity: windows into population history and trait architecture. Nat extraordinary diversity of immune genes. bioRxiv. 2018;347344. https://doi. Rev Genet. 2018;19(4):220–34. https://doi.org/10.1038/nrg.2017.109. org/10.1101/347344 . 43. Ceballos FC, Hazelhurst S, Ramsay M. Assessing runs of homozygosity: a 64. Miller YI, Choi S-H, Wiesner P, Bae YS. The SYK side of TLR4: Signalling comparison of SNP Array and whole genome sequence low coverage data. mechanisms in response to LPS and minimally oxidised LDL. Br J BMC Genomics. 2018;19(1):106. https://doi.org/10.1186/s12864-018-4489-0. Pharmacol. 2012;167(5):990–9. https://doi.org/10.1111/j.1476-5381.2012.02 44. Luigi-Sierra MG, Cardoso TF, Martínez A, Pons A, Bermejo LA, Jordana J, 097.x. et al. Low genome-wide homozygosity in 11 Spanish ovine breeds. Anim 65. Mócsai A, Ruland J, Tybulewicz VLJ. The SYK tyrosine kinase: a crucial player Genet. 2019;50(5):501–11. https://doi.org/10.1111/age.12832. in diverse biological functions. Nat Rev Immunol. 2010;10(6):387–402. 45. Comeron JM, Ratnappan R, Bailin S. The many landscapes of recombination https://doi.org/10.1038/nri2765. in Drosophila melanogaster. PLoS Genet. 2012;8(10):e1002905. https://doi. 66. Hou X, Lin L, Xing W, Yang Y, Duan X, Li Q, et al. Spleen tyrosine kinase org/10.1371/journal.pgen.1002905. regulates mammary epithelial cell proliferation in mammary glands of dairy 46. Smukowski CS, Noor MAF. Recombination rate variation in closely related cows. J Dairy Sci. 2016;99(5):3858–68. https://doi.org/10.3168/jds.2015-10118. species. Heredity. 2011;107(6):496–508. https://doi.org/10.1038/hdy.2011.44. 67. Wang X, Xu S, Gao X, Ren H, Chen J. Genetic polymorphism of TLR4 gene 47. Ohgami RS, Campagna DR, McDonald A, Fleming MD. The Steap proteins and correlation with mastitis in cattle. J Genet Genomics. 2007;34(5):406–12. are metalloreductases. Blood. 2006;108(4):1388–94. https://doi.org/10.1182/ https://doi.org/10.1016/S1673-8527(07)60044-7. blood-2006-02-003681. 68. Todd RF, Petty HR. β2(CD11/CD18) integrins can serve as signaling partners 48. Gomes IM, Maia CJ, Santos CR. STEAP proteins: from structure to for other leukocyte receptors. J Lab Clin Med. 1997;129(5):492–8. https://doi. applications in cancer therapy. Mol Cancer Res. 2012;10(5):537–87. https:// org/10.1016/S0022-2143(97)90003-2. doi.org/10.1158/1541-7786.MCR-11-0281. 69. Sokol CL, Luster AD. The chemokine system in innate immunity. Cold 49. Simmons RM, Erikson DW, Kim J, Burghardt RC, Bazer FW, Johnson GA, et al. Spring Harb Perspect Biol. 2015;7:1–20. Insulin-like growth factor binding protein-1 in the ruminant uterus: potential 70. Brenaut P, Lefèvre L, Rau A, Laloë D, Pisoni G, Moroni P, et al. Contribution endometrial marker and regulator of conceptus elongation. Endocrinology. of mammary epithelial cells to the immune response during early stages of 2009;150(9):4295–305. https://doi.org/10.1210/en.2009-0060. a bacterial infection to Staphylococcus aureus. Vet Res. 2014;45(1):16. https:// 50. Hammon HM, Zbinden Y, Sauerwein H, Breier BH, Blum JW, Donkin SS. The doi.org/10.1186/1297-9716-45-16. response of the hepatic insulin-like growth factor system to growth 71. Marsland BJ, Bättig P, Bauer M, Ruedl C, Lässing U, Beerli RR, et al. CCL19 hormone and dexamethasone in calves. J Endocrinol. 2003;179(3):427–35. and CCL21 induce a potent proinflammatory differentiation program in https://doi.org/10.1677/joe.0.1790427. licensed dendritic cells. Immunity. 2005;22(4):493–505. https://doi.org/10.101 51. Rupp R, Clément V, Piacere A, Robert-Granié C, Manfredi E. Genetic 6/j.immuni.2005.02.010. parameters for milk somatic cell score and relationship with production and 72. Yoshida H, Hunter CA. The immunobiology of interleukin-27. Annu Rev udder type traits in dairy Alpine and Saanen primiparous goats. J Dairy Sci. Immunol. 2015;33(1):417–43. https://doi.org/10.1146/annurev-immunol- 2011;7(7):3629–34. https://doi.org/10.3168/jds.2010-3694. 032414-112134. 52. Paape MJ, Poutrel B, Contreras A, Marco JC, Capuco AV. Milk somatic cells 73. Petes C, Odoardi N, Plater SM, Martin NL, Gee K. IL-27 amplifies cytokine and lactation in small ruminants. J Dairy Sci. 2001;84:E237–44. https://doi. responses to Gram-negative bacterial products and Salmonella typhimurium org/10.3168/jds.S0022-0302(01)70223-8. infection. Sci Rep. 2018;8(1):13704. https://doi.org/10.1038/s41598-018-32 53. Lerondelle C, Richard Y, Issartial J. Factors affecting somatic cell counts in 007-y. goat milk. Small Rumin Res. 1992;8(1-2):129–39. https://doi.org/10.1016/ 74. Bai X, Zheng Z, Liu B, Ji X, Bai Y, Zhang W. Whole blood transcriptional 0921-4488(92)90014-U. profiling comparison between different milk yield of Chinese Holstein cows using RNA-seq data. BMC Genomics. 2016;17(S7):512. https://doi.org/10.11 54. Kirin M, McQuillan R, Franklin CS, Campbell H, Mckeigue PM, Wilson JF. 86/s12864-016-2901-1. Genomic runs of homozygosity record population history and 75. Li Q, Liang R, Li Y, Gao Y, Li Q, Sun D, et al. Identification of candidate consanguinity. PLoS ONE. 2010;5(11):e13996. https://doi.org/10.1371/journal. genes for milk production traits by RNA sequencing on bovine liver at pone.0013996. different lactation stages. BMC Genet. 2020;21(1):72. https://doi.org/10.1186/ 55. Silió L, Rodríguez MC, Fernández A, Barragán C, Benítez R, Óvilo C, et al. s12863-020-00882-y. Measuring inbreeding and inbreeding depression on pig growth from 76. Schennink A, Bovenhuis H, Léon-Kloosterziel KM, Van Arendonk JAM, Visker pedigree or SNP-derived metrics. J Anim Breed Genet. 2013;130(5):349–60. MHPW. Effect of polymorphisms in the FASN, OLR1, PPARGC1A, PRL and https://doi.org/10.1111/jbg.12031. STAT5A genes on bovine milk-fat composition. Anim Genet. 2009;40(6):909– 56. Pemberton TJ, Szpiech ZA. Relationship between deleterious variation, 16. https://doi.org/10.1111/j.1365-2052.2009.01940.x. genomic autozygosity, and disease risk: insights from the 1000 genomes project. Am J Hum Genet. 2018;102(4):658–75. https://doi.org/10.1016/j.a jhg.2018.02.013. 57. Rokouei M, Vaez Torshizi R, Moradi Shahrbabak M, Sargolzaei M, Sørensen AC. Monitoring inbreeding trends and inbreeding depression for economically important traits of Holstein cattle in Iran. J Dairy Sci. 2010; 93(7):3294–302. https://doi.org/10.3168/jds.2009-2748. 58. Miglior F, Burnside EB, Dekkers JCM. Nonadditive genetic effects and inbreeding depression for somatic cell counts of Holstein cattle. J Dairy Sci. 1995;78(5):1168–73. https://doi.org/10.3168/jds.S0022-0302(95)76734-0. 59. Menzies PI, Ramanoon SZ. Mastitis of sheep and goats. Vet Clin North Am Food Anim Pract. 2001;17(2):333–58. https://doi.org/10.1016/S0749-0720(1 5)30032-3. 60. Spielman D, Brook BW, Briscoe DA, Frankham R. Does inbreeding and loss of genetic diversity decrease disease resistance? Conserv Genet. 2004;5(4): 439–48. https://doi.org/10.1023/B:COGE.0000041030.76598.cd. 61. Reid JM, Arcese P, Keller LF, Elliott KH, Sampson L, Hasselquist D. Inbreeding effects on immune response in free-living song sparrows (Melospiza melodia). Proc Biol Sci. 2007;274(1610):697–706. https://doi.org/10.1098/ rspb.2006.0092.