Skip to main content
Download PDF

Genomic patterns of homozygosity and inbreeding depression in Murciano-Granadina goats

Journal of Animal Science and Biotechnology volume 13, Article number: 35 (2022) Cite this article

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.

Background

Inbreeding is defined as the mating of individuals that are related to each other more closely than the average relationship within the concerned population [1]. In livestock species, the magnitude of inbreeding has been traditionally measured through genealogical information [2]. However, pedigree-based estimates are affected by the depth of the pedigree [2] because founders are assumed to be unrelated and non-inbred [3]. Consequently, inbreeding produced by distant ancestors not included in the pedigree is systematically ignored [4]. Another disadvantage of quantifying inbreeding from pedigree data is that it provides bare expectations about the fraction of the genome which is identical-by-descent (IBD) [3]. With the advent of high-density arrays of single nucleotide polymorphisms (SNPs), it has become possible to estimate genomic inbreeding coefficients which circumvent these limitations [5]. Indeed, important advantages of genomic inbreeding coefficients over their genealogical counterparts are: (i) higher accuracy to differentiate among individuals within the same pedigree, since variation due to Mendelian sampling is captured [4], (ii) higher accuracy to quantify shared ancestry of genetic haplotypes [4], and (iii) the ability to map inbreeding to specific genomic regions [6]. Different types of genomic inbreeding coefficients have been implemented. While inbreeding coefficients based on the proportion of homozygous SNPs (FHOM) just reflect identity-by-state (IBS) allele-sharing proportions [6], coefficients (FROH) based on measuring the fraction of the genome covered by runs of homozygosity (ROH) estimate IBD allele sharing [4,5,6,7,8], making possible to disentangle recent from ancient inbreeding [3,4,5,6,7,8,9].

The increase of inbreeding might have adverse consequences on the fitness of livestock populations due to the loss of genetic variability, which can entail a long-term reduction of genetic variance (due to the fixation of alleles) and, consequently, a slowing down of the rate of response to selection in breeding schemes [10,11,12]. Moreover, incremented levels of inbreeding might reduce the mean phenotypic performance of livestock populations, a phenomenon known as inbreeding depression (reviewed by Leroy [2]). Although inbreeding depression is particularly intense for fitness and reproduction traits [11], there is broad evidence that it also decreases dairy and growth performances [2, 13,14,15]. Besides, susceptibility to certain diseases, such as mastitis, is increased in inbred animals [16, 17]. In Holstein cattle, a 1% increase of inbreeding is expected to cause a reduction of $22–24 of lifetime net income per individual [18], while in sheep the average economic loss per ewe amounts to $17 for moderate inbreeding and $36 when inbreeding is close to 50% [18].

Several studies have used genomic methods to determine the levels of inbreeding in goat populations with a broad geographic distribution [19,20,21]. A recent investigation carried out by Bertolini et al. [19] revealed that short ROH (< 3 Mb) are particularly abundant in worldwide goat populations. Moreover, five regions on caprine chromosomes (CHI) 11, 12, and 18 contained ROH hotspots that overlapped with signatures of selection [19]. The majority of goat breeds analysed by Bertolini et al. [19] displayed low levels of inbreeding (FROH < 0.10), with the only exception of certain local breeds with small population sizes (e.g. Dutch Landrace goats) as well as of breeds with insular origins (e.g. Icelandic and Malagasy goats) which happened to be highly inbred [19, 20]. Despite the fact that inbreeding depression can have adverse effects on the profitability of farmers and animal breeders, very few studies have investigated its consequences on goat production [22, 23]. In this regard, Marete et al. [22] and Deroide et al. [23] estimated, with genealogical methods, the effect of inbreeding depression on the production of Kenya Alpine and Murciano-Granadina goats, respectively, and they found that in both populations such effect was negligible.

The goals of the current work were: (i) to measure the levels of inbreeding in a Murciano-Granadina resource population by using different genomic coefficients, and (ii) to use this information to infer the impact of inbreeding depression on dairy phenotypes recorded in this population.

Methods

Animal material and phenotyping

The animal material comprised 1040 Murciano-Granadina female goats from 15 farms located in the autonomous region of Andalusia (Spain). Murciano-Ganadina is a local Spanish breed officially created in 1975 by the crossbreeding of Murciano and Granadina goats [24]. Currently, it has a census of 115,105 heads (2020), and its remarkable adaptability to harsh environments as well as its good milking performance (mean of 586 kg/lactation; 5.1% of fat and 3.6% of protein in milk) have made it a very popular breed in Spain and other countries (https://www.mapa.gob.es/es/ganaderia/ temas/zootecnia/razas-ganaderas/razas/catalogo-razas/caprino/murciano-granadina/).

Blood samples were extracted from goats using vacuum tubes coated with EDTA K3 anticoagulant and stored at − 20 °C until processing. Phenotypic records for milk yield and composition traits were recorded in the framework of the selection program of the Murciano-Granadina goat breed [24]. Only phenotypes corresponding to the first parity (recorded between the years 2009 and 2017) were taken into consideration. The following phenotypes were recorded in 817 goats: milk yield measured in kilograms at 210 days (MY210), 240 days (MY240) and 305 days (MY305), the natural logarithm of the somatic cell count divided by 1000 (lnSCC, to convert this value into a somatic cell count please use the formula: elnSCC × 103 cells/mL), fat percentage (FP), protein percentage (PP) and lactose percentage (LP). Milk composition traits were standardised to a lactation of 210 days. Summary statistics of phenotypic records are displayed in Table 1.

Table 1 Summary statistics of seven milk production and composition traits recorded in 817 Murciano-Granadina goats
Full size table

Generation of high throughput genotypic data

The isolation of genomic DNA was carried out following a salting-out protocol [25]. Three mL of whole blood were mixed with 4 volumes of Red Cell Lysis Solution (Tris-HCl 10 mmol/L, pH = 6.5; EDTA 2 mmol/L; Tween 20 1%), and this mixture was centrifuged at 2000 × g. The supernatant was discarded and the pellet was resuspended in 3 mL of lysis buffer (Tris-HCl 200 mmol/L, pH = 8, EDTA 30 mmol/L, SDS 1%; NaCl 250 mmol/L) plus 100 μL proteinase K (20 mg/mL) and incubated for 3 h at 55 °C. The lysate was chilled, and 1 mL of ammonium acetate 10 mol/L was added to it. After centrifugation at 2000 × g for 10 min, the supernatant (~ 4 mL) was transferred to a new tube with 3 mL of isopropanol 96%, and this mixture was centrifuged at 2000 × g for 3 min. The resulting DNA pellet was washed with 3 mL of ethanol 70% and an additional centrifugation step at 2000 × g for 1 min was performed. The DNA pellet was dried at room temperature, and it was subsequently resuspended in 1 mL of TE buffer (Tris-HCl 10 mmol/L, EDTA 1 mmol/L, pH = 8). Murciano-Granadina goats were genotyped with the Goat SNP50 BeadChip (Illumina Inc., San Diego, CA, USA) by following the instructions of the manufacturer. The goat ARS1 genome [26] was used as reference for inferring the genomic location of the SNPs, and the position and the name of each SNP were updated using the PLINK 1.9 software [27]. Only individuals with at least 95% of SNPs with genotype calls were taken into consideration. Moreover, only SNPs meeting the following requirements were used in the downstream analyses: (i) mapping to autosomes, (ii) displaying a minor allele frequency of 0.05 or higher, (iii) not deviating significantly (P < 0.00001) from the Hardy-Weinberg equilibrium, and (iv) with a genotype call rate over 98%. Data were filtered using PLINK 1.9 [27]. In addition, a principal component analysis (PCA) was carried out with PLINK 1.9 [27] in order to assess population structure. The visualisation of such results was based on the first two components of the PCA.

Detection of runs of homozygosity

The definition of ROH followed six criteria: (i) the minimum length of ROH is 1 Mb, (ii) a ROH must contain at least 15 SNPs, (iii) the density of SNPs per ROH was set to at least 1 SNP every 100 kb, (iv) the maximum distance between consecutive SNPs is 250 Kb, (v) one heterozygous position per ROH is allowed, and (vi) one missing position per ROH is allowed. Said criteria were established based on the density of the genotyping panel, with a mean distance between consecutive SNPs of 51.73 Kb and a mean number of 19.35 SNPs/Mb (Additional file 1: Table S1), and several of them are based on the ROH definition established by the AdaptMap Consortium [19, 20]. Runs of homozygosity were identified with the PLINK 1.9 software [27] using a sliding window of 50 SNPs.

Analysis of the genomic patterns of homozygosity

The proportion of homozygosity per site was estimated as the ratio between the number of animals with homozygous genotypes for a particular SNP divided by the number of animals genotyped for that SNP. A sliding window encompassing 30 SNPs was designed to estimate the average of this ratio, and chromosomal patterns of homozygosity were visualised as Manhattan plots using R [28].

Genomic coverage and distribution of ROH were also investigated. The patterns of ROH size and distribution along the genome were analysed and plotted using R [28]. Genomic regions in which ROH are prevalent, the so-called ROH hotspots, were identified by measuring the proportion of animals that harbour a particular SNP occurring within a ROH with regard to the total number of animals genotyped for that SNP. Genomic regions containing the top 1% SNPs most commonly associated with ROH were classified as ROH hotspots [29, 30]. Both highly homozygous regions and ROH hotspots were compared, and overlapping segments were identified. Taking as reference the ARS1 goat genome [26], genes mapping to these overlapping segments were extracted with the Biomart tool of Ensembl [31]. Database for Annotation, Visualization and Integration Discovery (DAVID) Bioinformatics Resources [32] and UniProt [33] were used in Gene ontology analysis to identify over-represented (enriched) gene ontology (GO) terms and KEGG pathways. Amongst other things, DAVID provides biological context for long lists of genes by assigning them to functionally related groups through the use of fuzzy clustering techniques [32]. To assess the significance of gene-enrichment in annotation terms, a Fisher Exact is employed. Only terms with Max. Prob ≤ 0.1, Min. Count ≥ 2 and P < 0.05 were considered as significantly enriched [32].

Estimation of inbreeding coefficients

Six inbreeding coefficients were estimated at the whole genome level based on genotypic data:

FHOM is defined as the proportion of SNPs with homozygous genotypes [6] and was estimated as:

$$ {F}_{HOM_i}=\frac{H_{O_i}}{S} $$

where HOi corresponds to the observed number of homozygous genotypes for all the SNPs for each individual i and S are the number of SNPs for which individual i has genotype data. It was computed from the output of the --het command of PLINK 1.9 software [27].

FROH was estimated as the proportion of the genome covered by ROH by using the following formula:

$$ {F}_{ROH_i}=\frac{\sum {L}_{ROH_i}}{L_{auto}} $$

where LROHi corresponds to the sum of the lengths of all ROH present in each individual i, and Lauto is the total length of the autosomal goat genome covered by SNPs [5]. The same mathematical expression was used to calculate the genomic coverage of ROH with sizes smaller (FROHShort) or larger (FROHLong) than 5 Mb. Such calculations were made to assess the relative importance of distant (FROHShort) versus recent (FROHLong) inbreeding [6], according to the size of the genome and SNP density (Table 1). Assuming that the length of the ROH segments follows an exponential distribution with a mean equal to \( \frac{1}{2}g \), where g corresponds to the number of generations to the closest common ancestor [34], and also assuming that goats have a recombination rate of approximately 1 cM/Mb [35], FROHShort indicates the inbreeding of an individual from 10 to 50 generations in the past, while FROHLong estimates the inbreeding from 1 to 9 generations in the past.

FIS coefficient of Wright [36] was calculated with the formula:

$$ {F}_{IS}=\frac{H_{E_i}-{H}_{O_i}}{H_{E_i}} $$

Where \( {H}_{E_i} \) and \( {H}_{O_i} \) are the expected and observed heterozygosities of the individual i. This coefficient was estimated from the output of the --hardy command of PLINK 1.9 [27].

FYANG coefficient was estimated from the diagonal of the matrix of genomic relationships of Yang based on the correlation between uniting gametes [37, 38]. It was calculated with the following formula:

$$ {F}_{YANG_i}=\frac{1}{S}\sum \limits_{k=1}^s\frac{x_k^2-\left(1+2{p}_k\right){x}_{k_i}+2{p}_k^2}{2{p}_k\left(1-{p}_k\right)} $$

where, xk is the genotype of the individual i for the SNP k, and pk is the frequency of the reference allele in the studied population. The command --ibc of PLINK 1.9 [27] was used to estimate it.

The coefficients FHOM and FROH were also estimated at the chromosomal level.

Inbreeding depression analyses

The effects of inbreeding depression on dairy traits were investigated using data from 817 goats with available phenotypes (MY210, MY240, MY305, lnSCC, FP, PP and LP). Analyses were performed with the REMLF90 software [39] to implement a restricted maximum likelihood (REML) analysis approach, in which the phenotypic values of each trait in each individual are regressed onto its inbreeding coefficient using a linear mixed model. These analyses were performed to quantify inbreeding depression at the whole genome scale as well as on a chromosome and regional basis. The model was fitted as follows:

$$ y= X\beta + Za+e $$

where y is the vector of observations for each phenotype, β is a vector of fixed effects, including farm (15 levels), year of birth (10 levels), and the linear regression on F as a covariate; a is the vector of additive genetic effects, e is the vector of random residual effects, and X and Z are incidence matrices relating fixed and random effects to observations.

The significance of the inbreeding effect on the analysed traits was determined by applying a two-tailed hypothesis test. A Z-statistic was estimated with the following general formula:

$$ Z=\frac{\overline{\mathbf{x}}-{\mu}_0}{s.e.} $$

Here, \( \overline{x} \) corresponds to the regression coefficient representing the effect of inbreeding over each trait, and μ0 is the coefficient of inbreeding corresponding to the null hypothesis (in this case is equal to 0), and s.e. is the standard error. The transformation of Z-scores into P-values was accomplished with the function pnorm() implemented on R [28].

The above analysis was performed by regressing each phenotype onto four genomic inbreeding coefficients (FHOM, FROH, FROHLong, and FROHShort). In order to detect genomic regions associated with inbreeding depression, analyses at the chromosomal level were performed for traits that in the whole-genome analysis were identified as significantly affected by inbreeding depression (P < 0.05). Following Saura et al. [6], inbreeding depression was finely mapped by dividing chromosomes into six segments and performing the analyses reported above in each segment. As previously explained, genes mapping to genomic regions associated with inbreeding depression for a specific trait were retrieved using Biomart [31], and their biological functions were assessed with UniProt [33] and David Bioinformatics Resources version 6.8 [32].

Results

Assessment of homozygosity patterns in Murciano-Granadina goats

A total of 46,689 SNPs and 1040 animals were selected to investigate the population structure and patterns of homozygosity of Murciano-Granadina goats. The PCA (Additional file 2: Fig. S1) indicated that goats clustered, to some extent, in accordance with their farm of origin. We detected 20,312 ROH that were classified as follows: 11,325 had sizes of 0–5 Mb, while 5470 (5–10 Mb), 2695 (10–20 Mb), 789 (20–50 Mb) and 33 (> 50 Mb) displayed sizes above 5 Mb. The mean number of ROH per category and per individual varied slightly across farms (Fig. 1A and B). The mean ROH number was 19.53 ± 11.89 per individual, with an average length of 6.15 ± 2.05 Mb. As depicted in Fig. 2, the majority of the individuals harboured less than 50 ROH, and ROH covered a small proportion of the genome (< 300 Mb, about 10% of the genome). Only 2% of the individuals showed a genomic ROH coverage > 500 Mb, and 1% harboured more than 50 ROH. Larger chromosomes encompassed a greater number of ROH when compared to the smaller ones (Additional file 3: Fig. S2), and the correlation coefficient between the number of ROH and chromosome length was 0.92 (P < 0.05).

Fig. 1
figure 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

Full size image
Fig. 2
figure 2

Relationship between ROH number and length in 1040 Murciano-Granadina goats genotyped with the Goat SNP50 BeadChip

Full size image

The genome-wide analysis of homozygosity, based on the proportion of homozygous individuals for each genotyped position, made it possible to detect 38 genomic regions with increased homozygosity (harbouring the top 1% of the most homozygous positions) that were scattered on 20 goat chromosomes (CHI), i.e. CHI 1–8, 11, 13–18, 20, 21, 24, 26 and 29 (Fig. 3A and Additional file 4: Table S2). Eight ROH hotspots mapping to CHI 1, 2, 4, 6, 14, 16, and 17 were identified (Fig. 3B and Additional file 5: Table S3). One region (i.e. CHI4:42,552,375–48,378,207 bp) was consistently detected in the genome-wide analysis of homozygosity and ROH. Sixty-six genes mapped to these regions (Additional file 6: Table S4) and a functional enrichment analysis revealed 15 GO terms significantly enriched at the nominal level (P < 0.05). Particularly significant were GO terms related with ferric and copper import into the cell (Additional file 7: Table S5).

Fig. 3
figure 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

Full size image

Estimation of inbreeding coefficients

Genomic inbreeding coefficients reached values of 0.601 ± 0.021 (FHOM); 0.053 ± 0.046 (FROH); 0.040 ± 0.041 (FROHLong); 0.014 ± 0.008 (FROHShort); − 0.016 ± 0.035 (FIS) and 0.023 ± 0.047 (FYANG). The magnitude and dispersion of these coefficients are shown in Fig. 4. The FHOM, FROH, FROHLong, FIS and FYANG coefficients were highly correlated, being especially high the correlations between FHOM and FROH, FROH and FROHLong (r = 0.99, P< 2.2 × 10−16) and between FHOM and FIS (r = − 1, P < 2.2 × 10−16). In contrast, FROHShort showed the weakest correlations with the remaining inbreeding coefficients (|r| = 0.33–0.64), although their statistical significance (P < 2.2 × 10−16) was very high (Table 2).

Fig. 4
figure 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

Full size image
Table 2 Pearson correlations between molecular inbreeding coefficients (F) estimated in 1040 Murciano-Granadina goats
Full size table

Measurement of inbreeding depression for dairy traits

The natural logarithmic transformation of the somatic cell count divided by 1000 (lnSCC) as well as milk yield at three different time points (MY210, MY240 and MY305) showed significant evidence of inbreeding depression when analysing FHOM and FROH. In contrast, no significant effects of inbreeding were identified for coefficients based either on short or long ROH. Increases of 0.1 units of FHOM and FROH coefficients involved lnSCC increments of 0.29 (P = 0.037) and 0.127 (P = 0.038) units, respectively (Table 3). At the chromosomal level, significant inbreeding depression for lnSCC was detected on CHI8 and CHI25 when regressed on either FHOM or FROH, while six additional chromosomes (CHI13, CHI14, CHI22, CHI24, CHI25 and CHI27) displayed inbreeding depression for this trait exclusively when it was regressed onto FHOM (Additional file 8: Table S6). Four regions containing 666 genes on chromosomes 8 (i.e. CHI8:37,557,623–56,336,435 bp and CHI8:75,115,244–93,894,055 bp) and 25 (i.e. CHI25: 82,419–7,143,084 bp and CHI25:21,429,255–28,572,340 bp) displayed a significant inbreeding depression for lnSCC based on FHOM (Additional file 9: Table S7 and Additional file 10: Table S8). After performing functional enrichment analysis, various GO terms and pathways involved in the immune response were significantly enriched only at the nominal level, including integrin domains (i.e. ITGAM, ITGAD, ITGAX, ITGAL) and genes with protein kinase activity chemotaxis and cell signalling activities (i.e. CCL19, CCL21, CCL24, CCL26, SYK and IL27). Several of these genes display functions related with innate immunity and inflammatory response e.g. lymphocyte chemotaxis, cellular response to interferon-γ, immunological synapse formation, positive regulation of chemotaxis, cellular response to interleukin-1, monocyte chemotaxis etc. (Additional file 11: Table S9).

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)
Full size table

With regard to the milk yield traits, an increment of 0.1 units of FHOM or FROH involved a decrease of 48.50 kg (P = 0.011) and 20.49 kg (P = 0.014) of milk for MY210; 51.02 kg (P = 0.017) and 21.30 kg (P = 0.023) for MY240; and 55.72 kg (P = 0.032) and 23.20 kg (P = 0.040) for MY305, respectively (Table 3). The analysis at the chromosomal level indicated significant inbreeding depression for MY210 and MY240 on CHI2, CHI3 and CHI11, while for MY305 inbreeding depression was significant only on CHI2 (Additional file 12: Table S10). The region on CHI2, comprised between 22,751,824 and 68,255,473 bp displayed significant inbreeding for all three milk yield traits (Additional file 13: Table S11). A total of 355 genes mapped to this genomic region (Additional file 14: Table S12). The functional enrichment analysis evidenced an overrepresentation of genes involved in mitochondrial and energetic processes, including genes from the PPAR signalling pathway i.e. CYP27A1, ACADL, DBI and ACSL9 (Additional file 15: Table S13).

Although there is no substantial overlap between regions associated with inbreeding depression and ROH hotspots, the CHI8:75–93.8 Mb region associated with inbreeding depression for lnSCC and the CHI2:45.5–68.25 Mb region associated with inbreeding depression for milk yield show positional concordance with two of the 38 regions displaying high homozygosity (i.e. CHI8:76.25–77.4 Mb and CHI2:56.13–57.17 Mb).

Discussion

Low inbreeding in Murciano-Granadina goats

For the majority of Murciano-Granadina goats, ROH number and total length were below 50 and 350 Mb, respectively. Moreover, short ROH (< 5 Mb) were more abundant than the medium or long ones (Fig. 1A). These patterns are pretty consistent with what has been observed by Bertolini et al. [19] in a worldwide sample of goat breeds. Indeed, Bertolini et al. [19] showed that goat breeds from Southern Europe had, on average, 49 ROH per individual while the genomic coverage per individual was 183.47 Mb. In contrast, breeds from Northern Europe showed higher levels of homozygosity with 98 ROH per individual and genomic coverage of 479.17 Mb.

The inbreeding coefficients FROH, FROHShort, and FROHLong of Murciano-Granadina goats were mainly in the range of 0 to 0.05. In their study, Bertolini et al. [19] reported that about 60% of a worldwide sample of goat breeds displayed low FROH coefficients (< 0.10), while the remaining ~ 30% and ~ 10% of breeds showed moderate (0.10 < FROH < 0.20) or high (> 0.20) FROH values. The patterns of low homozygosity that we have observed in Murciano-Granadina goats contrast strongly with what has been reported in certain local breeds, such as Mallorquina, Pyrenean, and Valdostana, which have undergone sharp population bottlenecks [19, 20]. Low inbreeding and homozygosity in the Murciano-Granadina breed are probably explained by its large census size (> 100,000 individuals registered in the herdbook), the absence of population bottlenecks, and its broad geographic distribution encompassing more than 4000 farms across Spain (https://www.mapa.gob.es/es/ganaderia/temas/zootecnia/razas-ganaderas/razas) and other countries. Noteworthy, the Murciano-Granadina breed was founded by crossing, during the 1970s, two Murciana and Granadina populations with different historical origins [40]. Although genetic differentiation between these two populations was weak [40], this admixture event probably contributed to increase the heterozygosity of the resulting Murciano-Granadina breed. Widespread use of artificial insemination in reproductive management and intensive selection were implemented in the Murciano-Granadina breed a few decades ago [24], so their impact on genetic diversity and inbreeding has probably been quite limited so far.

While correlations between FHOM, FROH, FROHLong, FIS and FYANG were high, FROHShort displayed the lowest correlations with the remaining inbreeding coefficients, in line with previous studies focused on other livestock species [6, 13, 41]. Short ROH are mainly derived from ancient inbreeding events [5, 42] and do not reflect the whole autozygosity of the sample. It is also possible that several of these homozygous tracks are identical by state and not by descent, being produced by a low recombination rate or high linkage disequilibrium in unrelated ancestors [30]. Besides, when working with medium density genotype arrays (e.g. 50 K SNPs) the detection of short ROH can become quite difficult [43], thus limiting the ability to infer the true proportion of short vs. long ROH in the genome.

Several ROH hotspots are detected in the genomes of Murciano-Granadina goats

We have identified several genomic regions in which ROH are particularly frequent (ROH hotspots). Similar patterns were found when the proportion of homozygous individuals was analysed at a genome-wide level i.e. 35.46% of the most common homozygous regions overlapped with ROH hotspots and these overlapping regions represented 35.42% of the total ROH hotspots. No positional coincidence was detected between ROH hotspots identified by us and those reported by Bertolini et al. [19] in a worldwide sample of goat populations. This finding agrees well with what has been previously observed in sheep [44]. Indeed, ROH hotspots are produced by factors such as positive selection and inversions suppressing recombination, that can differ substantially from population to population [42]. For instance, the characterisation of the products of 5860 female meioses in Drosophila melanogaster by genotyping more than 100 million SNPs made it possible to detect 106,964 recombination events displaying a remarkable degree of intra-specific variation [45]. Factors such as GC content, gene density, distribution of simple repeats and transposable elements, structural variation, and the presence of diverse poorly-characterised sequence motifs might explain the regional variation of the recombination rate across individuals and populations [46].

ROH hotspots contain genes with diverse functions

Regarding the gene content of genomic segments co-localising with both high homozygosity regions and ROH hotspots (CHI4:42,552,375–48,378,207 bp), the functional enrichment analysis highlighted several gene ontology terms with nominally significant enrichment (P-value < 0.05) (Additional file 7: Table S5). From this list of genes, it is worth emphasising STEAP1, STEAP2 and STEAP4 that encode homonymous metalloreductases which facilitate the cellular uptake of iron and copper [47]. These proteins modulate the effects of intracellular oxidative stress and inflammation and are involved in multiple biological pathways related with molecular trafficking in the endocytic and exocytic pathways, metabolism, control of cell proliferation and apoptosis and tumour progression [48]. We also found genes with metabolic functions such as insulin-like growth factor binding protein 3 (IGFBP3) and insulin-like growth factor binding protein 1 (IGFBP1) that participate in the growth and postnatal development of cattle [49, 50]. These results are quite concordant with those reported by Mastrangelo and colleagues [30], who showed that ROH islands identified in Italian bovine breeds contained genes with heterogeneous functions related to milk production, reproduction, immune response, and resistance/susceptibility to infection and diseases.

Effect of inbreeding depression on milk traits

The lnSCC in the Murciano-Granadina population under study averaged 6.25 ± 0.93 units, which is higher than values reported in primiparous goats from the Alpine (5.09 ± 1.36 units) and Saanen (5.32 ± 1.19 units) breeds [51]. Compared with cows and ewes, goats display higher numbers of somatic cells in milk. Indeed, the apocrine nature of milk secretion in goats increases the proportion of cytoplasmatic particles in milk, a feature, that depending on the measurement method of choice,  might increase the somatic cell count [52]. Besides, somatic cell count is modulated by many factors including the occurrence of bacterial infections, stress, oestrous cycle phase, diet etc. [53]. According to our results (Table 3), inbreeding depression increased lnSCC, a feature that is considered adverse because high lnSCC values are often associated with subclinical and clinical mastitis [51, 52].

The magnitude of inbreeding depression for lnSCC estimated from FROHShort strongly differed from estimates based on the other inbreeding coefficients (Table 3), a finding consistent with the moderate (|r| = 0.33–0.64) correlations between FROHShort and other molecular inbreeding coefficients (Table 2). Noteworthy, recent rather than ancestral inbreeding is the main cause of inbreeding depression in mammalian populations [54, 55]. Besides, long stretches of homozygosity usually contain a higher proportion of deleterious mutations than the shorter ones because they are more recent, so deleterious variation has not been yet purged by purifying selection [56].

In Murciano-Granadina goats, a previous study performed by Deroide et al. [23] reported a low percentage of inbreeding (average F = 0.24%). Milk production showed a positive quadratic correlation with inbreeding levels, but such effect was not significant. Deroide et al. [23] also reported that milk fat and dry extract contents experienced a slight increase due to inbreeding. In our study, the dairy trait mostly affected by inbreeding depression was lnSCC. Consistently, Doekes et al. [13] reported that a 1% increase of FROH involved a 0.86 ± 0.28 unit increase in somatic cell score (days 150 through to 400) recorded in Dutch Holstein-Friesian dairy cattle. In Iranian cattle, individuals with high inbreeding coefficients tended to have higher somatic cell scores than animals with low inbreeding coefficients [57], and similar results have been reported for Canadian Holstein cattle [58]. Doekes et al. [13] indicated that ancient inbreeding was the main contributor to inbreeding depression for somatic cell score, although such effect was not significant. Somatic cell score is an indicator of the health status of the mammary gland and substantial increases are observed in individuals suffering from mastitis [59]. Inbreeding has been reported to significantly reduce resistance against pathogens in multiple organisms [60,61,62], so the significant inbreeding depression observed for lnSCC in Murciano-Granadina goats might be explained, at least in part, by the weakening of the immune defences of the mammary gland. Thus, homozygosity for deleterious mutations might result in the partial or total inactivation of genes related with immunity, and low variability might also compromise the effectiveness of the immune response [63].

Genomic regions associated with inbreeding depression for lnSCC contain several genes related with immunity

When we investigated which enriched clusters are present in the set of genes mapping to chromosomes 8 (37–56 Mb and 75–93 Mb) and 25 (0.082–7 Mb and 21–28 Mb) regions associated with inbreeding depression for lnSCC (Additional file 10: Table S8), we found several genes assigned to gene ontologies highly connected with immunity, e.g. integrin-mediated signalling pathway, lymphocyte chemotaxis, monocyte chemotaxis, immunological synapse formation, chemokine activity, etc. (Additional file 11: Table S9). The spleen tyrosine kinase (SYK) protein forms part of the integrins cluster and maps to CHI8: 86,755,291–86,861,895 bp (Additional file 11: Table S9). One of the functions of the SYK molecule is to stimulate the phosphorylation of Toll-like receptor 4 (TLR4) [64], which recognises bacterial lipopolysaccharide and induces inflammatory and immune responses [65]. This gene has been described as highly variable in cattle [66], and many TLR4 polymorphisms and haplotypes have been associated with milk somatic cell count and susceptibility to mastitis [67]. Moreover, the SYK protein induces the recognition of pathogens and cell adhesion and platelet activation [65], and it also affects the proliferation of mammary epithelial cells at several stages of the milking cycle [66]. In the same enriched gene ontologies, we have detected the integrin subunit αM gene (ITGAM also known as CD11b) which maps to CHI25: 27,221,164–27,264,350 bp and encodes a receptor for lipopolysaccharide [68]. Signalling mediated by TLR4 activates the synthesis of ITGAM/CD11b, which is essential for the migration and adhesion of polymorphonuclear leukocytes to infection sites [68].

As previously indicated, genes related to chemotaxis were significantly enriched at the nominal level (P-value < 0.05) (Additional file 11: Table S9). This functional category is mainly represented by chemokine genes, such as chemokine ligand 19 (CCL19), 21 (CCL21), 24 (CCL24), 26 (CCL26) and 27 (CCL27). Chemokines are essential for the development of the innate immune response since they orchestrate and control the migration of the immune cells (macrophages, monocytes, neutrophils, etc.) to sites of infection [69]. During the first stages of mastitis, chemokines contribute to the stimulation of the cellular immune response against the invading pathogen until acute-phase proteins are expressed [70]. Marsland et al. [71] described how chemokines CCL19 and CCL21 participate in the maturation of dendritic cells, allowing them to leverage the T cell response. These chemokines also enhance the migration of leukocytes through lymph and blood circulation and stimulate the production of pathogen-induced proinflammatory 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 with both pro and anti-inflammatory effects, thus enhancing the immune response and, at the same time, preventing tissue damage caused by inflammation [72]. Infections caused by Gram-negative bacteria induce IL27 production, and this cytokine interacts with monocytes increasing TLR4 expression and enhancing the LPS-induced inflammatory response [73]. Moreover, IL27 has an autocrine effect on macrophages and monocytes resulting in the amplification of the inflammatory response via cytokine secretion [72].

The genomic region displaying inbreeding depression for milk yield (MY210, MY240 and MY305) was significantly enriched at the nominal level (P-value < 0.05) with genes associated with multiple unrelated biological processes (Additional file 15: Table S13). The PPAR signalling pathway (P-value = 0.032) influences milk production and composition in cattle. Bai and collaborators [74] described an overrepresentation of genes from the PPAR signalling pathway in cows with high daily milk yield in comparison with low yielders. Besides, genes from this pathway are upregulated in cows in the lactation peak when compared to those in the dry period [75], suggesting a role of these genes not only in the determinism of fat composition [76] but also of milk yield.

Conclusions

Murciano-Granadina goats display low levels of inbreeding (mean FROH = 0.053 ± 0.046), a finding consistent with the large census size and demographic history of this breed. Four genomic regions on CHI8:37,557,623–56,336,435, CHI8:75,115,244–93,894,055 bp, CHI25:82,419–7,143,084 and CHI25:21,429,255–28,572,340 bp were associated with inbreeding depression for lnSCC. Moreover, one region on CHI2:22,751,824–68,255,473 was consistently associated with inbreeding depression for three milk yield traits (MY210, MY240 and MY310). Genes encoding integrins, chemokines and pathogen recognition receptors, which play relevant roles in the elicitation of innate immune responses against microbes, mapped to regions associated with inbreeding depression for lnSCC. These results suggest that keeping inbreeding to a minimum, through an adequate reproductive management, might be a useful approach to decrease the incidence of mastitis in Murciano-Granadina goats.

Availability of data and materials

Goat SNP50 BeadChip genotypes and milk production phenotypes are accessible at Figshare (https://doi.org/10.6084/m9.figshare.18095825).

Abbreviations

IBD:

Identical-by-descent

SNPs:

Single nucleotide polymorphisms

IBS:

Identity-by-state

ROH:

Runs of homozygosity

CHI:

Chromosomes

FP:

Fat percentage

PP:

Protein percentage

LP:

Lactose percentage

PCA:

Principal component analysis

DAVID:

Annotation Visualization and Integration Discovery

GO:

Gene ontology

REML:

Restricted maximum likelihood

SYK:

Spleen tyrosine kinase

References

  1. Lush JL. The consequences and measurement of inbreeding. In: Animal breeding plans. Ames: The Iowa State College Press; 1943. 259–288.

  2. Leroy G. Inbreeding depression in livestock species: review and meta-analysis. Anim Genet. 2014;45(5):618–28. https://doi.org/10.1111/age.12178.

    Article  CAS  PubMed  Google Scholar 

  3. Kardos M, Luikart G, Allendorf FW. Measuring individual inbreeding in the age of genomics: marker-based measures are better than pedigrees. Heredity. 2015;115(1):63–72. https://doi.org/10.1038/hdy.2015.17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Keller MC, Visscher PM, Goddard ME. Quantification of inbreeding due to distant ancestors and its detection using dense single nucleotide polymorphism data. Genetics. 2011;189(1):237–49. https://doi.org/10.1534/genetics.111.130922.

    Article  PubMed  PubMed Central  Google Scholar 

  5. McQuillan R, Leutenegger AL, Abdel-Rahman R, Franklin CS, Pericic M, Barac-Lauc L, et al. Runs of homozygosity in European populations. Am J Hum Genet. 2008;83(3):359–72. https://doi.org/10.1016/j.ajhg.2008.08.007.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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 using genome-wide data. Genet Sel Evol. 2015;47(1):1. https://doi.org/10.1186/s12711-014-0081-5.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Broman KW, Weber JL. Long homozygous chromosomal segments in reference families from the Centre d’Etude du Polymorphisme Humain. Am J Hum Genet. 1999;65(6):1493–500. https://doi.org/10.1086/302661.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Pemberton TJ, Absher D, Feldman MW, Myers RM, Rosenberg NA, Li JZ. Genomic patterns of homozygosity in worldwide human populations. Am J Hum Genet. 2012;91(2):275–92. https://doi.org/10.1016/j.ajhg.2012.06.014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Purfield DC, Berry DP, McParland S, Bradley DG. Runs of homozygosity and population history in cattle. BMC Genet. 2012;13(1):70. https://doi.org/10.1186/1471-2156-13-70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Kristensen TN, Sørensen AC. Inbreeding - lessons from animal breeding, evolutionary biology and conservation genetics. Anim Sci. 2005;80(2):121–33. https://doi.org/10.1079/ASC41960121.

    Article  Google Scholar 

  11. Charlesworth D, Willis JH. The genetics of inbreeding depression. Nat Rev Genet. 2009;10(11):783–96. https://doi.org/10.1038/nrg2664.

    Article  CAS  PubMed  Google Scholar 

  12. Howard JT, Pryce JE, Baes C, Maltecca C. Invited review: inbreeding in the genomics era: inbreeding, inbreeding depression, and management of genomic variability. J Dairy Sci. 2017;100(8):6009–24. https://doi.org/10.3168/jds.2017-12787.

    Article  CAS  PubMed  Google Scholar 

  13. Doekes HP, Veerkamp RF, Bijma P, De Jong G, Hiemstra SJ, Windig JJ. Inbreeding depression due to recent and ancient inbreeding in Dutch Holstein-Friesian dairy cattle. Genet Sel Evol. 2019;52(1):64. https://doi.org/10.1186/s12711-020-00583-1.

    Article  CAS  Google Scholar 

  14. Pryce JE, Haile-Mariam M, Goddard ME, Hayes BJ. Identification of genomic regions associated with inbreeding depression in Holstein and Jersey dairy cattle. Genet Sel Evol. 2014;46(1):71. https://doi.org/10.1186/s12711-014-0071-7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Thompson JR, Everett RW, Hammerschmidt NL. Effects of inbreeding on production and survival in Holsteins. J Dairy Sci. 2000;83(9):2131–8. https://doi.org/10.3168/jds.S0022-0302(00)75096-X.

    Article  CAS  PubMed  Google Scholar 

  16. Sørensen AC, Madsen P, Sørensen MK, Berg P. Udder health shows inbreeding depression in Danish Holsteins. J Dairy Sci. 2006;89(10):4077–82. https://doi.org/10.3168/jds.S0022-0302(06)72452-3.

    Article  PubMed  Google Scholar 

  17. Croquet C, Mayeres P, Gillon A, Vanderick S, Gengler N. Inbreeding depression for global and partial economic indexes, production, type, and functional traits. J Dairy Sci. 2006;89(6):2257–67. https://doi.org/10.3168/jds.S0022-0302(06)72297-4.

    Article  CAS  PubMed  Google Scholar 

  18. Ercanbrack SK, Knight AD. Effects of inbreeding on reproduction and wool production of Rambouillet, Targhee, and Columbia ewes. J Anim Sci. 1991;69(12):4734–44. https://doi.org/10.2527/1991.69124734x.

    Article  CAS  PubMed  Google Scholar 

  19. Bertolini F, Cardoso TF, Marras G, Nicolazzi EL, Rothschild MF, Amills M. Genome-wide patterns of homozygosity provide clues about the population history and adaptation of goats. Genet Sel Evol. 2018;50(1):59. https://doi.org/10.1186/s12711-018-0424-8.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Cardoso TF, Amills M, Bertolini F, Rothschild M, Marras G, Boink G, et al. Patterns of homozygosity in insular and continental goat breeds. Genet Sel Evol. 2018;50(1):56. https://doi.org/10.1186/s12711-018-0425-7.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Nandolo W, Mészáros G, Banda LJ, Gondwe TN, Lamuno D, Mulindwa HA, et al. Timing and extent of inbreeding in African goats. Front Genet. 2019;10:537. https://doi.org/10.3389/fgene.2019.00537.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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.

    Google Scholar 

  23. Deroide CAS, Jacopini LA, Delgado JV, Léon JM, Brasil LHA, Ribeiro MN. Inbreeding depression and environmental effect on milk traits of the Murciano-Granadina goat breed. Small Rumin Res. 2016;134:44–8. https://doi.org/10.1016/j.smallrumres.2015.12.008.

    Article  Google Scholar 

  24. Delgado JV, Landi V, Barba CJ, Fernández J, Gómez MM, Camacho ME, et al. Murciano-Granadina goat: a Spanish local breed ready for the challenges of the twenty-first century. In: Simões J, Gutiérrez C, editors. Sustainable goat production in adverse environments: volume II. Springer: Cham; 2018. p. 205–19. https://doi.org/10.1007/978-3-319-71294-9_15.

    Chapter  Google Scholar 

  25. Miller SA, Dykes DD, Polesky HF. A simple salting out procedure for extracting DNA from human nucleated cells. Nucleic Acids Res. 1988;16(3):1215. https://doi.org/10.1093/nar/16.3.1215.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Bickhart DM, Rosen BD, Koren S, Sayre BL, Hastie AR, Chan S, et al. Single-molecule sequencing and chromatin conformation capture enable de novo reference assembly of the domestic goat genome. Nat Genet. 2017;49(4):643–50. https://doi.org/10.1038/ng.3802.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Chang CC, Chow CC, Tellier LCAM, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience. 2015;4(1):7. https://doi.org/10.1186/s13742-015-0047-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. ISBN 3–900051–07-0, URL http://www.R-project.org/

  29. Purfield DC, McParland S, Wall E, Berry DP. The distribution of runs of homozygosity and selection signatures in six commercial meat sheep breeds. PLoS ONE. 2017;12(5):e0176780. https://doi.org/10.1371/journal.pone.0176780.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Mastrangelo S, Tolone M, Sardina MT, Sottile G, Sutera AM, Di Gerlando R, et al. Genome-wide scan for runs of homozygosity identifies potential candidate genes associated with local adaptation in Valle del Belice sheep. Genet Sel Evol. 2017;49(1):84. https://doi.org/10.1186/s12711-017-0360-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Kinsella RJ, Kähäri A, Haider S, Zamora J, Proctor G, Spudich G, et al. Ensembl BioMarts: A hub for data retrieval across taxonomic space. Database. 2011;2011:bar030.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57. https://doi.org/10.1038/nprot.2008.211.

    Article  CAS  Google Scholar 

  33. Bateman A. UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 2019;47(D1):D506–15. https://doi.org/10.1093/nar/gky1049.

    Article  CAS  Google Scholar 

  34. Thompson EA. Identity by descent: variation in meiosis, across genomes, and in populations. Genetics. 2013;194(2):301–26. https://doi.org/10.1534/genetics.112.148825.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Dumont BL, Payseur BA. Evolution of the genomic rate of recombination in mammals. Evolution. 2008;62(2):276–94. https://doi.org/10.1111/j.1558-5646.2007.00278.x.

    Article  CAS  PubMed  Google Scholar 

  36. Wright S. The genetical structure of populations. Ann Eugenics. 1951;15(4):323–54. https://doi.org/10.1111/j.1469-1809.1949.tb02451.x.

    Article  CAS  Google Scholar 

  37. Villanueva B, Fernández A, Saura M, Caballero A, Fernández J, Morales-González E, et al. The value of genomic relationship matrices to estimate levels of inbreeding. Genet Sel Evol. 2021;53(1):42. https://doi.org/10.1186/s12711-021-00635-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Yang J, Benyamin B, McEvoy BP, Gordon S, Henders AK, Nyholt DR, et al. Common SNPs explain a large proportion of the heritability for human height. Nat Genet. 2010;42(7):565–9. https://doi.org/10.1038/ng.608.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Misztal I, Tsuruta S, Strabel T, Auvray B, Druet T, Lee DH. BLUPF90 and related programs (BGF90). In Proceedings of the 7th World Congress on Genetics Applied to Livestock Production: 19-23 August 2002; Montpellier. 2002:28.07.

  40. Martínez A, Rocha L, Quiroz J, Delgado BJ. Estudio de la diversidad genética intrarracial de la cabra Murciano-Granadina con microsatélites de ADN. Arch Zootec. 2007;56:417–20.

    Google Scholar 

  41. Biscarini F, Mastrangelo S, Catillo G, Senczuk G, Ciampolini R. Insights into genetic diversity, runs of homozygosity and heterozygosity-rich regions in Maremmana semi-feral cattle using pedigree and genomic data. Animals. 2020;10(12):2285. https://doi.org/10.3390/ani10122285.

    Article  PubMed Central  Google Scholar 

  42. Ceballos FC, Joshi PK, Clark DW, Ramsay M, Wilson JF. Runs of homozygosity: windows into population history and trait architecture. Nat Rev Genet. 2018;19(4):220–34. https://doi.org/10.1038/nrg.2017.109.

    Article  CAS  PubMed  Google Scholar 

  43. Ceballos FC, Hazelhurst S, Ramsay M. Assessing runs of homozygosity: a comparison of SNP Array and whole genome sequence low coverage data. BMC Genomics. 2018;19(1):106. https://doi.org/10.1186/s12864-018-4489-0.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Luigi-Sierra MG, Cardoso TF, Martínez A, Pons A, Bermejo LA, Jordana J, et al. Low genome-wide homozygosity in 11 Spanish ovine breeds. Anim Genet. 2019;50(5):501–11. https://doi.org/10.1111/age.12832.

    Article  CAS  PubMed  Google Scholar 

  45. Comeron JM, Ratnappan R, Bailin S. The many landscapes of recombination in Drosophila melanogaster. PLoS Genet. 2012;8(10):e1002905. https://doi.org/10.1371/journal.pgen.1002905.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Smukowski CS, Noor MAF. Recombination rate variation in closely related species. Heredity. 2011;107(6):496–508. https://doi.org/10.1038/hdy.2011.44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Ohgami RS, Campagna DR, McDonald A, Fleming MD. The Steap proteins are metalloreductases. Blood. 2006;108(4):1388–94. https://doi.org/10.1182/blood-2006-02-003681.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Gomes IM, Maia CJ, Santos CR. STEAP proteins: from structure to applications in cancer therapy. Mol Cancer Res. 2012;10(5):537–87. https://doi.org/10.1158/1541-7786.MCR-11-0281.

    Article  CAS  Google Scholar 

  49. Simmons RM, Erikson DW, Kim J, Burghardt RC, Bazer FW, Johnson GA, et al. Insulin-like growth factor binding protein-1 in the ruminant uterus: potential endometrial marker and regulator of conceptus elongation. Endocrinology. 2009;150(9):4295–305. https://doi.org/10.1210/en.2009-0060.

    Article  CAS  PubMed  Google Scholar 

  50. Hammon HM, Zbinden Y, Sauerwein H, Breier BH, Blum JW, Donkin SS. The response of the hepatic insulin-like growth factor system to growth hormone and dexamethasone in calves. J Endocrinol. 2003;179(3):427–35. https://doi.org/10.1677/joe.0.1790427.

    Article  CAS  PubMed  Google Scholar 

  51. Rupp R, Clément V, Piacere A, Robert-Granié C, Manfredi E. Genetic parameters for milk somatic cell score and relationship with production and udder type traits in dairy Alpine and Saanen primiparous goats. J Dairy Sci. 2011;7(7):3629–34. https://doi.org/10.3168/jds.2010-3694.

    Article  CAS  Google Scholar 

  52. Paape MJ, Poutrel B, Contreras A, Marco JC, Capuco AV. Milk somatic cells and lactation in small ruminants. J Dairy Sci. 2001;84:E237–44. https://doi.org/10.3168/jds.S0022-0302(01)70223-8.

    Article  Google Scholar 

  53. Lerondelle C, Richard Y, Issartial J. Factors affecting somatic cell counts in goat milk. Small Rumin Res. 1992;8(1-2):129–39. https://doi.org/10.1016/0921-4488(92)90014-U.

    Article  Google Scholar 

  54. Kirin M, McQuillan R, Franklin CS, Campbell H, Mckeigue PM, Wilson JF. Genomic runs of homozygosity record population history and consanguinity. PLoS ONE. 2010;5(11):e13996. https://doi.org/10.1371/journal.pone.0013996.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Silió L, Rodríguez MC, Fernández A, Barragán C, Benítez R, Óvilo C, et al. Measuring inbreeding and inbreeding depression on pig growth from pedigree or SNP-derived metrics. J Anim Breed Genet. 2013;130(5):349–60. https://doi.org/10.1111/jbg.12031.

    Article  CAS  PubMed  Google Scholar 

  56. Pemberton TJ, Szpiech ZA. Relationship between deleterious variation, 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.ajhg.2018.02.013.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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.

    Article  CAS  PubMed  Google Scholar 

  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.

    Article  CAS  PubMed  Google Scholar 

  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(15)30032-3.

    Article  CAS  PubMed  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  PubMed  Google Scholar 

  62. Smallbone W, van Oosterhout C, Cable J. The effects of inbreeding on disease susceptibility: Gyrodactylus turnbulli infection of guppies, Poecilia reticulata. Exp Parasitol. 2016;167:32–7. https://doi.org/10.1016/j.exppara.2016.04.018.

  63. Siljestam M, Rueffler C. Heterozygote advantage can explain the extraordinary diversity of immune genes. bioRxiv. 2018;347344. https://doi.org/10.1101/347344 .

  64. Miller YI, Choi S-H, Wiesner P, Bae YS. The SYK side of TLR4: Signalling mechanisms in response to LPS and minimally oxidised LDL. Br J Pharmacol. 2012;167(5):990–9. https://doi.org/10.1111/j.1476-5381.2012.02097.x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Mócsai A, Ruland J, Tybulewicz VLJ. The SYK tyrosine kinase: a crucial player in diverse biological functions. Nat Rev Immunol. 2010;10(6):387–402. https://doi.org/10.1038/nri2765.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Hou X, Lin L, Xing W, Yang Y, Duan X, Li Q, et al. Spleen tyrosine kinase regulates mammary epithelial cell proliferation in mammary glands of dairy cows. J Dairy Sci. 2016;99(5):3858–68. https://doi.org/10.3168/jds.2015-10118.

    Article  CAS  PubMed  Google Scholar 

  67. Wang X, Xu S, Gao X, Ren H, Chen J. Genetic polymorphism of TLR4 gene and correlation with mastitis in cattle. J Genet Genomics. 2007;34(5):406–12. https://doi.org/10.1016/S1673-8527(07)60044-7.

    Article  CAS  PubMed  Google Scholar 

  68. Todd RF, Petty HR. β2(CD11/CD18) integrins can serve as signaling partners for other leukocyte receptors. J Lab Clin Med. 1997;129(5):492–8. https://doi.org/10.1016/S0022-2143(97)90003-2.

    Article  CAS  PubMed  Google Scholar 

  69. Sokol CL, Luster AD. The chemokine system in innate immunity. Cold Spring Harb Perspect Biol. 2015;7:1–20.

    Article  Google Scholar 

  70. Brenaut P, Lefèvre L, Rau A, Laloë D, Pisoni G, Moroni P, et al. Contribution of mammary epithelial cells to the immune response during early stages of a bacterial infection to Staphylococcus aureus. Vet Res. 2014;45(1):16. https://doi.org/10.1186/1297-9716-45-16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Marsland BJ, Bättig P, Bauer M, Ruedl C, Lässing U, Beerli RR, et al. CCL19 and CCL21 induce a potent proinflammatory differentiation program in licensed dendritic cells. Immunity. 2005;22(4):493–505. https://doi.org/10.1016/j.immuni.2005.02.010.

    Article  CAS  PubMed  Google Scholar 

  72. Yoshida H, Hunter CA. The immunobiology of interleukin-27. Annu Rev Immunol. 2015;33(1):417–43. https://doi.org/10.1146/annurev-immunol-032414-112134.

    Article  CAS  PubMed  Google Scholar 

  73. Petes C, Odoardi N, Plater SM, Martin NL, Gee K. IL-27 amplifies cytokine responses to Gram-negative bacterial products and Salmonella typhimurium infection. Sci Rep. 2018;8(1):13704. https://doi.org/10.1038/s41598-018-32007-y.

  74. Bai X, Zheng Z, Liu B, Ji X, Bai Y, Zhang W. Whole blood transcriptional profiling comparison between different milk yield of Chinese Holstein cows using RNA-seq data. BMC Genomics. 2016;17(S7):512. https://doi.org/10.1186/s12864-016-2901-1.

    Article  PubMed  PubMed Central  Google Scholar 

  75. Li Q, Liang R, Li Y, Gao Y, Li Q, Sun D, et al. Identification of candidate genes for milk production traits by RNA sequencing on bovine liver at different lactation stages. BMC Genet. 2020;21(1):72. https://doi.org/10.1186/s12863-020-00882-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Schennink A, Bovenhuis H, Léon-Kloosterziel KM, Van Arendonk JAM, Visker MHPW. Effect of polymorphisms in the FASN, OLR1, PPARGC1A, PRL and STAT5A genes on bovine milk-fat composition. Anim Genet. 2009;40(6):909–16. https://doi.org/10.1111/j.1365-2052.2009.01940.x.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Many thanks to CAPRIGRAN for carrying out phenotype recording and blood sample collection in Murciano-Granadina goats.

Funding

This research was funded by the European Regional Development Fund (FEDER)/Ministerio de Ciencia e Innovación - Agencia Estatal de Investigación/Project Reference grant: PID2019-105805RB-I00 and by the CERCA Programme/Generalitat de Catalunya. We also acknowledge the support of the Spanish Ministerio de Ciencia e Innovación for the Center of Excellence Severo Ochoa 2020–2023 (CEX2019–000902-S) grant awarded to the Centre for Research in Agricultural Genomics (CRAG, Bellaterra, Spain). We also acknowledge the support of the CERCA programme of the Generalitat de Catalunya. Dailu Guan was funded by a PhD fellowship from the China Scholarship Council (CSC). Maria Luigi-Sierra was funded with a PhD fellowship Formación de Personal Investigador (BES-C-2017-079709) awarded by the Spanish Ministry of Economy and Competitivity.

Author information

Authors and Affiliations

  1. Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Universitat Autònoma de Barcelona, 08193, Bellaterra, Spain

    María Gracia Luigi-Sierra, Dailu Guan & Marcel Amills

  2. Departamento de Mejora Genética Animal, INIA, Carretera de la Coruña km 7,5, 28040, Madrid, Spain

    Almudena Fernández & María Saura

  3. Departamento de Genética, Universidad de Córdoba, 14071, Córdoba, Spain

    Amparo Martínez, Juan Vicente Delgado & Javier Fernández Álvarez

  4. Department of Veterinary Medicine, University of Bari ‘‘Aldo Moro”, 62 per Casamassima km. 3, 70010, Valenzano, SP, Italy

    Vincenzo Landi

  5. Group of Research in Ruminants (G2R), Department of Animal and Food Science, Universitat Autònoma de Barcelona (UAB), Bellaterra, Barcelona, Spain

    Francesc Xavier Such

  6. Departament de Ciència Animal i dels Aliments, Facultat de Veterinària, Universitat Autònoma de Barcelona, 08193, Bellaterra, Spain

    Jordi Jordana & Marcel Amills

Authors
  1. María Gracia Luigi-Sierra
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Almudena Fernández
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Amparo Martínez
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Dailu Guan
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Juan Vicente Delgado
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Javier Fernández Álvarez
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Vincenzo Landi
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Francesc Xavier Such
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Jordi Jordana
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. María Saura
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Marcel Amills
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

MA, JJ, JVD and VL designed the study. JFA, JVD and AM coordinated all 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 the cooperation of MS and AF. All authors read and approved the content of the paper.

Corresponding authors

Correspondence to María Saura or Marcel Amills.

Ethics declarations

Ethics approval and consent to participate

Blood collection is a routine procedure performed by CAPRIGRAN, so it does not require approval by the Ethics Committee on Animal and Human Experimentation of the Universitat Autònoma de Barcelona.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Supplementary Information

Additional file 1: Table S1.

Number and density of SNPs per chromosome in a population of 1040 Murciano-Granadina goats genotyped with the Goat SNP50 BeadChip.

Additional file 2: Figure S1.

Principal component analysis of 1040 Murciano-Granadina goats distributed in 15 farms (each farm is indicated with a different colour).

Additional file 3: Figure S2.

Number of ROH per chromosome (represented as yellow bars, left axis) and the percentage of each chromosome covered by ROH (represented by a red line, right axis) in 1040 Murciano-Granadina goats.

Additional file 4: Table S2.

Genomic regions associated with the top 1% of homozygosity for each SNP marker in a population of 1040 Murciano-Granadina goats genotyped with the Goat SNP50 BeadChip.

Additional file 5: Table S3.

ROH hotspots identified in the genomes of 1040 Murciano-Granadina goats.

Additional file 6: Table S4.

Genes mapping to regions consistently identified as ROH hotspots and regions with high homozygosity.

Additional file 7: Table S5.

Functional enrichment analysis of genes mapping to regions consistently identified as ROH hotspots and regions with high homozygosity.

Additional file 8: Table S6.

Inbreeding depression estimates (s.e: standard error) per goat chromosome (CHI) for the natural logarithm of the somatic cell count divided by 1000 (lnSCC) and their 95% confidence intervals (C.I).

Additional file 9: Table S7.

Inbreeding depression estimates (s.e: standard error) for the natural logarithm of the somatic cell count divided by 1000 (lnSCC) in specific regions of goat chromosomes (CHI) 8 and 25 and their 95% confidence intervals (C.I).

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) regions associated with inbreeding depression for the natural logarithm of the somatic cell count divided by 1000 (lnSCC).

Additional file 11: Table S9.

Functional enrichment analysis of genes 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 natural logarithm of the somatic cell count divided by 1000 (lnSCC).

Additional file 12: Table S10.

Inbreeding depression estimates (s.e: standard error) per goat chromosome (CHI) for milk yield at 210, 240 and 305 days (MY210, MY240 and MY305, measured in kg) and their 95% confidence intervals (C.I).

Additional file 13: Table S11.

Inbreeding depression estimates (s.e: standard error) for milk yield at 210, 240 and 305 days (MY210, MY240 and MY305, measured in kg) in specific regions of goat chromosomes (CHI) 2 and their 95% confidence intervals (C.I).

Additional file 14: Table S12.

Genes mapping to goat chromosome (CHI) 2: 22–68 Mb region associated with inbreeding depression for milk yield at 210, 240 and 305 days (MY210, MY240 and MY305).

Additional file 15: Table S13.

Functional enrichment analysis of genes from the goat chromosome (CHI) 2: 22–68 Mb region associated with inbreeding depression for milk yield at 210, 240 and 305 days (MY210, MY240 and MY305).

Rights and permissions

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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Luigi-Sierra, M.G., Fernández, A., Martínez, A. et al. Genomic patterns of homozygosity and inbreeding depression in Murciano-Granadina goats. J Animal Sci Biotechnol 13, 35 (2022). https://doi.org/10.1186/s40104-022-00684-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40104-022-00684-5

Keywords

Journal of Animal Science and Biotechnology

ISSN: 2049-1891