Estimation of genetic parameters and the effect of inbreeding on dairy traits in Istrian sheep

Abstract

Pedigree analysis is an important tool for assessing genetic variability and detecting heritage patterns in domestic animals. This study aimed to (1) assess genetic variability, (2) estimate genetic parameters for daily milk yield (DMY), fat content (FC), and protein content (PC), and (3) estimate the impact of inbreeding on the aforementioned dairy traits in Istrian sheep breed. The generation interval was 3.8 years. The average family size per sire and dam was 13.5 and 1.6, respectively. There were 923 female and 282 male founders (f), 105 male and 137 female effective founders (f_e), and 72 male and 95 female effective ancestors (f_a). A total of 28 sires and 34 dams explained 50% of the genetic variability. The effective number of male and female founder genomes (f_ne) was 33 and 53, respectively. Low f_e/f (0.37 and 0.15), high f_e/f_a (1.45 and 1.44), and low f_ne/f_e (0.45 and 0.39) for males and females indicated overuse of a small number of animals in reproduction and substantial loss of founder alleles. The heritabilities of MY, FC, and PC were 0.09, 0.03, and 0.06, respectively. Genetic correlations between DMY and FC, DMY and PC, and FC and PC were −0.22, −0.27, and 0.70, respectively. Regression of phenotypes on inbreeding coefficient ($F$) within the animal genetic model revealed negligible inbreeding depression for all traits (p<.01; a 1% increase in $F$ decreased DMY, FC, and PC by 0.002 kg/day, 0.01%, and 0.002%, respectively). Regression of breeding values on $F$ was insignificant for all traits (p>.05).

HIGHLIGHTS

• Genetic diversity, genetic parameters and the impact of inbreeding on dairy traits were estimated using pedigree information.

• Results indicate recent loss of genetic variability, low heritability and negligible impact of inbreeding on dairy traits.

• Genomic optimum contribution selection is the best solution to achieve genetic gain and slow down the loss of genetic variability in this population.

Keywords:

• Istrian sheep
• dairy traits
• genetic diversity
• genetic parameters
• inbreeding

Introduction

Dairy sheep farming has a long tradition in Mediterranean countries. Sheep cheese and curd are the main dairy products, usually priced above those of other milk. The Istrian sheep originates from the Istrian peninsula and has been bred for centuries in semi-arid climatic conditions with long, dry, hot summers and mild winters. Nowadays, the breed is also widespread in the Karst region of neighbouring countries Slovenia and Italy (Piasentier et al. 2002). Originally, the breed was classified as a dual-purpose breed (milk and meat), but more recently, with average milk production of 160.9 kg (CAAF 2020) and selection towards dairy traits, it is classified more like a dairy breed (Mioč et al. 2012). The milk recording system started in early 2000 following ICAR guidelines and the breeding programme was officially established in 2011 (Mioč et al. 2011). BLUP genetic evaluation for dairy traits has been conducted for about a decade now.

Strong selection pressure usually results in a substantial genetic gain, but intensive use of a limited number of genetically superior animals leads to higher inbreeding rates and loss of genetic variability (Falconer and Mckay 1996; Mrode 2014). The impact of selection on genetic variability can be assessed from the pedigree by calculating rates of inbreeding, co-ancestry, and effective population size (de Rochambeau et al. 2000; Goyache et al. 2003), as well as the contribution of founders, ancestors, effective founders, and founder genomes (Maignel et al. 1996; Boichard et al. 1997). Estimates of some of the above parameters have been reported in numerous studies (e.g. Eteqadi et al. 2014; Rodriguez-Ramilo et al. 2019; Barbosa et al. 2020; Domínguez-Viveros et al. 2020; Hashemi and Hossein-Zadeh 2020).

Genetic parameters such as heritability and genetic correlation represent important information to predict response to selection (Carta et al. 2009). Test‒day records from routine milking controls in modern sheep breeding programs usually represent the phenotypes for both, the estimation of genetic parameters and prediction of breeding values (BVs). Genetic parameters for dairy traits have been reported prevalently for French and Spanish dairy sheep populations such as Lacune (Barillet and Boichard 1994), Corsican, Red-Faced Manech, Black-Faced Manech, and Basco-Béarnaise breeds (Astruc et al. 2002); Churra (Baro et al. 1994; El-Saied et al. 1998), Latxa and Manchega sheep (Serrano et al. 2001; Legarra et al. 2003). However, there are also reports for breeds from other countries such as Slovenia: Bovec, Improved Bovec, and Istrian Pramenka (Brežnik 1999; Komprej et al. 2009); Czechia: East Friesian, Lacaune, Sumavska, and crossbreeds (Bauer et al. 2012); and Slovakia: Tsigai and Improved Valachian (Oravcová et al. 2005; Oravcova 2014), and Lacune (Oravcova 2007).

When kinship is neglected, selection and mating based solely on BLUP results have a tendency to increase the rate of inbreeding, which can be accompanied by inbreeding depression and deleterious recessive alleles in the homozygous state (Falconer and Mckay 1996). Despite numerous reports on the impact of inbreeding on growth traits in sheep and dairy traits in cattle, there are no reports on the impact of inbreeding on dairy traits in sheep. However, Murphy et al. (2017) hypothesised that due to the determined negative correlation between inbreeding and retained heterosis, inbreeding depression on dairy traits may also be present in sheep as well. In many cases, the non-significant effects of inbreeding on performance traits in sheep may have been the result of low pedigree completeness and a certain number of animals without phenotypic records. To bridge this gap in their studies, some authors decided to evaluate the effect of inbreeding on BVs (Barbosa et al. 2020).

The impact of non-genetic factors on dairy traits has been extensively studied in Istrian sheep so far (Vrdoljak et al. 2012; Kasap et al. 2019, 2021), but without addressing the issue of population genetic parameters. The objectives of this study were to (1) assess genetic diversity, (2) estimate genetic parameters using the multi-trait repeatability test–day model for daily milk yield (DMY), fat content (FC), and protein content (PC), and (3) estimate the effect of inbreeding on the abovementioned dairy traits.

Material and methods

Phenotypic and pedigree data

Phenotypic records for dairy traits (DMY, FC, and PC) and pedigree information were provided by the Ministry of Agriculture. Phenotypic records were collected according to alternate AT4 schemes (morning/evening system) following the ICAR guidelines (ICAR 2018). Initially, there were 51,405 test–day records collected in the period from February 2006 to August 2020. Before the final inferential statistical analysis, the original data set was ‘pruned’ as described by Spehar et al. (2020). Summary statistics for dairy traits used in the inferential statistical analysis are presented in Table 1.

Table 1. Descriptive statistics for dairy traits.

Trait	N	Mean	SD	Min	Max
Daily milk yield (kg)	37,703	1.015	0.523	0.204	2.938
Fat content (%)	37,607	7.39	1.76	2.00	14.99
Protein content (%)	37,613	5.97	0.72	3.00	9.00

The pedigree file was created using all available relationships to phenotyped ewes (5,071). After pruning of the initial file (7,050 animals), a total of 6,853 animals were included in the pedigree (Table 2). The proportion of base animals (both ancestors unknown) was ∼14%. Of the ∼86% of non–base animals, ∼90% had both parents known, 3% only sire known, and 6% only dam known. Pedigree analysis was used to characterise the population structure based on the following parameters:

Table 2. Pedigree structure.

Item	N	Proportion (%)
Ewes with records	5,071
Maximum number of generations traced back	10
Non–base animals	5,885	85.9
– both parents known	5,354
– only sire known	205
– only dam known	326
Base animals	968	14.1
Total number of animals	6,853

Demographic structure

1. Generation interval (GI) ‒ defined as the average age of the parents at the birth of their selected offspring (Falconer and Mckay 1996) and was calculated considering the four pathways: sire–son, sire–daughter, dam–son and dam–daughter.

2. Family size ‒ the number of offspring of an individual that becomes breeding individuals in the next generation (Falconer and Mckay 1996). The genetic contribution of the parents is usually not the same, leading to variation in family size and consequently an increase in the rate of inbreeding.

Pedigree quality

1. Number of fully traced generations (NTG) ‒ represents the number of generations separating an individual from its furthest ancestor.

2. Maximum number of complete generations (NCG) ‒ identifies the furthest generation with two known ancestors.

3. Number of equivalent complete generations (NECG) ‒ expresses the sum of all known ancestors counting how many generations have been traced. It is calculated for those individuals having at least one parent known following Maignel et al. (1996): (1) $$\mathit{NECG}=\sum _{i=1}^{n_j}\frac{1}{2^{\mathrm{g}}ij}\mathrm{ }$$(1)

   where $n_j\mathrm{ }$ is the number of ancestors $i$ of individual $j,$ and $g_{ij\mathrm{ }}$ is the number of generations between individual $j$ and its ancestor $i.$ Since it is the sum of the proportion of known ancestors over all generations, one-half is added for each known parent, one-fourth is added for each known grandparent, and so on.

4. Pedigree completeness ‒ is related to the number of generations that separate the individual from its most distant ancestor and the number of equivalent full generations and was assessed by calculating the pedigree completeness index (PCI) based on five generations as proposed by MacCluer et al. (1983). This parameter summarises the proportion of known ancestors in each ascending generation as follows: (2) $${\mathit{PCI}}_{\mathrm{animal}}=\frac{{2C}_{\mathrm{sire}}+{2C}_{\mathrm{dam}}}{C_{\mathrm{sire}}+C_{\mathrm{dam}}}\mathrm{ }$$(2)

   where $C_{\mathrm{sire}}\mathrm{ }$ and $C_{\mathrm{dam}}$ are contributions of sire and dam ancestors calculated using the following equation: (3) $$C=\frac{1}{d}\sum _{i=1}^d{a}_i\mathrm{ }$$(3)

   where $a_i$ is the proportion of known ancestors in generation $i$ and $d$ is the number of generations considered.

In order to provide enough information for accurate estimation of population genetic parameters, only animals with sufficient pedigree information (born from 2010 to 2018 with PCI ≥ 0.5 and NECG ≥ 3) were retained in the reference population (Figure 1).

Figure 1. Average pedigree completeness index (PCI) by year of birth for animals in the pedigree (the green area corresponds to the reference population).

Parameters based on probability of gene origin

1. Number of founders (f) ‒ founder is an animal with unknown genetic connections to other animals in the pedigree except for its own progeny (Lacy 1989).

2. Effective number of founders ($f_e$) ‒ because some founders have been used more intensely and therefore contribute more to the current population than other founders, the number of equally contributing founders represents founders that would be expected to produce the same genetic diversity as assessed in the population (Lacy 1989; Boichard et al. 1997). Essentially,$\mathrm{ }{f}_e\mathrm{ }$ denotes the number of animals with an equal genetic contribution, i.e. those that have contributed a similar amount of genetic material to the population (Pedrosa et al. 2010). The genetic contributions of the founders are independent and sum to one. It can be calculated as follows: (4) $$f_e={\left[\sum _{i=1}^{N_f}{q}_i^2\right]}^{-1}\mathrm{ }$$(4)

   where $q_i$ is the genetic contribution of the $i_{th}$ founder to the reference population. When founders contribute unequally, the effective number of founders is smaller than the actual number. The ratio $f_e/f$ indicates the degree of unequal founder contribution and this ratio is generally below 1.

3. Effective number of ancestors ($f_a)$ ‒ the minimum number of animals (not necessarily founders), explaining the complete genetic diversity of a population (Boichard et al. 1997). It accounts for the losses of genetic variability caused by the unbalanced use of reproductive individuals leading to bottlenecks. The $f_a$ is lower than the $f_e$ because it accounts for bottlenecks in the pedigree and is also less sensitive to missing pedigrees (Sørensen et al. 2005). The $f_a\mathrm{ }$ is calculated as follows: (5) $$f_a={\left[\sum _{i=1}^a{p}_i^2\right]}^{-1}\mathrm{ }$$(5)

   where $p_i$ is the marginal genetic contribution of ancestor $i.$ For the calculation of the $f_a\mathrm{ }$ one must consider only the marginal contribution of an ancestor. The ancestors contributing the most were chosen in an iterative process, where at round $n,$ the $p_i,$ is defined as the genetic contribution of ancestor $i,$ not yet explained by the $n-1$ already selected ancestors (Boichard et al. 1997). Ancestors have a large marginal contribution to the reference population if their genes have passed through many descendants, e.g. a sire of sons with many sons selected. We used ancestors explaining 50% of the genetic diversity in the case when the upper and lower limits did not show any difference for two decimal places. If the $f_a$/$f_e$ is close to unity, the population has been stable in terms of number of effectively contributing animals. If the $f_e$ is larger than the $f_a,$ bottlenecks have played a role in population development (Sørensen et al. 2005).

4. Effective number of founder genomes ($f_{ne}$) ‒ considers unequal contributions of founders, bottlenecks, and random loss of alleles due to genetic drift (Lacy 1989). It is calculated as: (6) $$f_{ne}={\left[\sum _{i=1}^{N_f}\frac{q_i^2}{r_i}\right]}^{-1}\mathrm{ }$$(6)

where $q_i\mathrm{ }$ is expected genetic contribution of the $\mathit{ith}$ founder and $r_i\mathrm{ }$ is expected proportion of founders $i\mathrm{\text{'}}s$ alleles which remain in the current population. The $f_{ne}$ is an indication of the random loss of alleles. Since alleles are lost in every generation, the effective number of founder genomes decreases in every generation and is therefore sensitive to depth of pedigree (Sørensen et al. 2005). The ratio $f_{ne}$/$f_e$ indicates the impact of drift where lower values indicates more loss of genetic diversity through drift.

Parameters based on probability of gene origin

1. Coefficient of inbreeding ($F$) ‒ it was determined by calculating $F$ using Cholesky factorisation of the relationship matrix as proposed by Meuwissen and Luo (1992). The inbreeding rate ($\Delta F$) was estimated by regressing $F$ on the generation interval (sliding window approach).

2. Coefficient of co-ancestry ‒ the $F$ of an individual is equal to half of the additive genetic relationship (AGR) between its parents or the coefficient of co-ancestry (Falconer and Mckay 1996; Mrode 2014). The co-ancestry of animals in one generation forecasts the average $F$ in the following generation (Falconer and Mckay 1996).

PEDIG (Boichard 2002) and optiSel package (Wellmann 2019) in the R programming environment (R Core Team 2020) were used to analyse the above parameters.

Covariance components estimation and prediction of breeding values

Covariance components were estimated by REML using analytical gradients in the VCE–6 (Groeneveld et al. 2008) based on multi-trait repeatability test–day model. The same model was used to predict BVs. The model included the following fixed class effects: parity (1–6), litter size (singles, twins), the season of lambing (58 classes), and flock (92 classes). Days in milk (DIM), age at lambing, and the $F$ were fitted as covariates. DIM was modelled using the Ali–Schaeffer lactation curve (Ali and Schaeffer 1987) nested within parity and litter size for DMY and within parity for FC and PC. Age at lambing was modelled as linear regression nested within parity in the model for DMY. The random part of the model was the same for all traits and included flock–test–day, permanent environmental effect within lactations, direct additive genetic effect, and residual.

The matrix notation of the model is as follows: (7) $$\mathrm{y}=\mathbf{Xb}+{\mathbf{Z}}_{\mathrm{c}}\mathbf{c}+{\mathbf{Z}}_{\mathrm{p}}\mathbf{p}+{\mathbf{Z}}_{\mathrm{a}}\mathbf{a}+\mathbf{e}\mathrm{ }$$(7) where: y is a vector of the phenotypic observations, X is an incidence matrix for the fixed effects; ${\mathbf{Z}}_{\mathrm{c}},$ ${\mathbf{Z}}_{\mathrm{p}},$ and ${\mathbf{Z}}_{\mathrm{a}},$ are incidence matrices for the flock–test–day, permanent environmental, and additive genetic effects, respectively; $\mathbf{b}\mathrm{ }$ is a vector of unknown fixed effects; $\mathbf{c},$ $\mathbf{p},$ and $\mathbf{a},$ are vectors of unknown random effects; $\mathbf{e}$ is a vector of residuals. The variance structure of random effects is: $\mathrm{Var}\mathrm{ }\left[\mathbf{c}\right]={\mathbf{I}}_{\mathrm{c}}{\mathrm{\sigma }}_{\mathrm{c}}^2,$ $\mathrm{Var}\mathrm{ }\left[\mathbf{p}\right]={\mathbf{I}}_{\mathrm{p}}{\mathrm{\sigma }}_{\mathrm{p}}^2,$ $\mathrm{Var}\mathrm{ }\left[\mathbf{e}\right]={\mathbf{I}}_{\mathrm{e}}{\mathrm{\sigma }}_{\mathrm{e}}^2,$ where $\mathbf{I}$ is an identity matrix for individual random effects and residual, $\mathit{Var}\mathrm{ }\left[\mathit{a}\right]=\mathit{A}{\sigma }_a^2,$ where $\mathit{A}\mathrm{ }$ is an additive genetic relationship matrix, and $\mathrm{Var}\mathrm{ }\left[\mathbf{y}\right]={\mathbf{Z}}_{\mathbf{a}}\mathbf{A}{\mathbf{Z}}_{\mathbf{a}}^{\mathrm{\text{'}}}{\mathrm{\sigma }}_{\mathrm{a}}^2+{\mathbf{Z}}_{\mathbf{c}}{\mathbf{Z}}_{\mathbf{c}}^{\mathrm{\text{'}}}{\mathrm{\sigma }}_{\mathrm{c}}^2+{\mathbf{Z}}_{\mathbf{p}}{\mathbf{Z}}_{\mathbf{p}}^{\mathrm{\text{'}}}{\mathrm{\sigma }}_{\mathrm{p}}^2+{\mathbf{I}}_{\mathbf{e}}{\mathrm{\sigma }}_{\mathrm{e}}^2.$ The assumption was that residuals, as well as trivial random effects, were not correlated.

Impact of inbreeding on dairy traits

The effect of inbreeding on phenotypes was obtained by regressing test day records on $F$ in the fixed part of the animal genetic model. The impact of inbreeding on dairy traits was additionally estimated by regressing BVs on $F.$

Results

The demographic structure of Istrian sheep is presented in Table 3. The average GI was 3.8 years. The GI differed between sires and dams. The longest GI was determined in the dam–daughter pathway which lasted 4.3 years, followed by dam–son (4.0 years). Higher GI for dam lines are often observed as females are usually kept in breeding for a longer period. The GI for sire–son and sire–daughter pathways were shorter and lasted 3.4 and 3.5 years, respectively. The maximum number of breeding progeny per sire and dam was 93 and 7, respectively. The average family size of sires and dams was 13.5 and 1.6, respectively.

Table 3. Demographic structure.

Item	N	Mean
Generation interval
Males for sire–son	333	3.4
Females for sire–daughter	2,408	3.5
Males dam–son	335	4.0
Females for dam–daughter	2,483	4.3
Family size	Maximum
Per sire	93	13.5
Per dam	7	1.6

Pedigree completeness of the reference population is presented in Table 4. NECG, NTG, NCG, and PCI were 4.12, 2.69, 6.74, and 0.88, respectively. Figure 1 shows PCI trends for the total and reference population. The steep increase in pedigree completeness observed from 2000 to 2010 was not surprising, as official pedigree recording in this population began in the early 2000s. The PCI for the reference population was above 80% throughout the period studied and exceeded 90% in the last two years. The correlation between the NECG and PCI was 0.96, confirming that they assess more or less the same thing.

Table 4. Pedigree quality of the reference population.

Item	Mean	SD	Min	Max
NECG	4.12	0.72	3	6.63
NTG	2.69	0.79	1	5
NCG	6.74	1.52	3	12
PCI	0.88	0.09	0.54	1.00

NECG: number of equivalent complete generations; NTG: number of fully traced generations; NCG: maximum number of complete generations; PCI: pedigree completeness index.

Estimates of the population parameters involved in assessing genetic variability are presented in Table 5. The number of female and male founders ($f$) was 923 and 282, and the number of effective founders ($f_e$) 105 and 137, respectively. The $f_a,$ which represents the minimum number of ancestors to allow a full explanation of the genetic diversity of the population, was 72 for males and 95 for females. A relatively small number of ancestors (28 males and 34 females) explained 50% of the genetic variability. The estimated $f_{ne}$ was 33 for males and 53 for females, indicating the presence of lower proportions of founder genes in the reference population. The $f_e/f$ ratio, which estimates genetic loss due to unequal contribution of the founders was low, especially for ewes (0.15), evidencing substantial historical overuse of some animals. The $f_e$/$f_a\mathrm{ }$ ratio was equally high for both genders (∼1.45) indicating a substantial loss of variability due to genetic bottleneck. The $f_{ne}$/$f_e$ ratio was similarly low for males (0.45) and females (0.39), additionally proving a substantial loss of founder alleles between generations.

Table 5. Pedigree-based estimates derived from the probabilities of gene origin.

Parameter	Male	Female
Number of founders ($f$)	282	923
Effective number of founders ($f_e$)	105	137
Effective number of ancestors ($f_a$)	72	95
N. of ancestors explaining 50% of genetic variability	28	34
Effective number of founder genomes ($f_{ne}$)	33	53
$f_e /f$ ratio, $f_e\mathrm{ }/\mathrm{ }{f}_a\mathrm{ }\mathrm{ratio},\mathrm{ }{f}_{ne}/\mathrm{ }{f}_e$ ratio	0.37/1.45/0.45	0.15/1.44/0.39

The proportion of inbred and non-inbred animals was 72.05% and 27.95%, respectively (Table 6). The highest proportion of animals (about 30%) had $F$<6.25%, while $F$ between 6.25% and 18.75% was found in about 30% of the animals. A high level of inbreeding ($F$>25%) was found in 10% of the population. The average $F$ in the reference population was 7.83%. The average coefficient of co-ancestry was 0.0116. The average rate of inbreeding per generation was 0.31% (Figure 2).

Figure 2. The number of all animals, number of inbred animals and average $F$ per generation.

Table 6. Distribution of inbreeding coefficient ($F$) by classes of inbreeding and descriptive statistics of $F$ in the reference population.

Coefficient of inbreeding ($F$)	Animals (N)	Animals (%)	Average
0	472	27.95
0.001<F<6.25	497	29.43
6.25≤F<12.5	256	15.16
12.5≤F<18.75	222	13.14
18.75≤F<25	77	4.56
F>25	165	9.77
Inbreeding (%)			7.83
Inbreeding of inbred animals (%)			10.87
Rate of inbreeding per year (%)			0.32
Coefficient of co-ancestry			0.0116

Heritabilities of MY (0.09 ± 0.004), FC (0.03 ± 0.002), and PC (0.06 ± 0.002) obtained from the multi-trait repeatability test–day model were low compared to commonly reported values for these traits in dairy sheep populations (Table 7). Additive genetic correlations between DMY and FC, and DMY and PC were −0.22 and −0.27, respectively (Table 8). The estimated additive genetic correlation between FC and PC was 0.70.

Table 7. Ratios of variance components in total variance (± standard errors).

Item	Daily milk yield (kg)	Fat content (%)	Protein content (%)
$c^2$	0.40 ± 0.003	0.35 ± 0.003	0.21 ± 0.002
$p^2$	0.19 ± 0.004	0.05 ± 0.003	0.08 ± 0.003
$h^2$	0.09 ± 0.004	0.03 ± 0.002	0.06 ± 0.002
$e^2$	0.32 ± 0.003	0.58 ± 0.004	0.66 ± 0.004

$c^2$ – flock-test day effect; $p^2$ – permanent environmental effect; $h^2$ – heritability; $e^2$ – residual.

Table 8. Genetic correlations between dairy traits.

Item	FC	PC
Flock-test-day
DMY	−0.13 ± 0.01	−0.07 ± 0.01
FC		0.20 ± 0.01
Permanent environmental effect
DMY	−0.14 ± 0.02	−0.20 ± 0.01
FC		0.54 ± 0.02
Additive genetic effect
DMY	−0.22 ± 0.02	−0.27 ± 0.02
FC		0.70 ± 0.02
Residual
DMY	−0.05 ± 0.01	−0.07 ± 0.01
FC		0.37 ± 0.004

*DMY – daily milk yield; FC – fat content; PC – protein content.

Regression of phenotypic values on the $F$ revealed negligible but significant inbreeding depression (p < .01) for all traits (Table 9). It was estimated that a 1% increase in $F$ decreased DMY, FC, and PC for 0.002 kg/day, 0.01%, and 0.002%, respectively. The regression of breeding values on $F$ was both negligible and insignificant (p > .05) for all traits analysed.

Table 9. Estimation of regression slope of $F$ on dairy traits.

Trait	Phenotype			Breeding values
	Estimate	Std Error	p	Estimate	Std Error	p
Daily milk yield (kg)	−0.0021	0.00038	<.0001	−0.0004	0.00147	>.05
Fat content (%)	−0.0079	0.00124	<.0001	−0.0005	0.00046	>.05
Protein content (%)	−0.0017	0.00051	<.0001	0.0014	0.00069	>.05

Discussion

Monitoring genetic variability in populations undergoing selection is an important step in livestock breeding programs and can be accomplished using a variety of ‘tools’. Nowadays, genomic data overtake the main role in this field (Meuwissen et al. 2020), however, pedigree-derived parameters still represent the major source of information in many studies and conservation breeding programs (Goyache et al. 2003; Gutiérrez et al. 2003). Estimation of genetic parameters related to diversity (effective number of founders, ancestors, inbreeding, etc.) and inheritance of traits (heritabilities and genetic correlations) from genealogical data significantly depends on the construction of the reference population. Only animals with sufficient ancestral information fairly contribute to the estimates of population-specific parameters and only they should compose the reference population (Boichard et al. 1997; Gutiérrez et al. 2003; Sørensen et al. 2005). Deep quality control of the overall pedigree in this population revealed a steep decline in the pedigree completeness with each ancestral generation which was not surprising since both, phenotype and pedigree recording, in this population started about two decades ago. However, the ‘poor’ pedigree quality parameters were related to the very old animals that were not included in further analysis. Animals born after 2010 showed a much better pedigree profile, especially after the additional pruning (NECG, NTG, NCG, and PCI were 4.12, 2.69, 6.74, and 0.88, respectively).

Even though it was not the main goal of the study, we estimated the generation interval (GI) because it represents an important population parameter for predicting response to selection and provides an indication about dynamics of generations change in the population under concern. Since increased selection intensity tends to reduce GI, it results in a more rapid loss of genetic variability (Boichard et al. 1997; Gutiérrez et al. 2003), it can also provide some insight about the level of genetic erosion in populations under selection. The average GI (3.8 years) fits well with estimates of this parameter in many dairy sheep populations (Barros et al. 2017; Hashemi and Hossein-Zadeh 2020) suggesting intermediate selection pressure. Taking into account selection for milk yield and the above-noted relationship between GI and loss of genetic variability, the shorter sire-progeny pathway leads us to the conclusion that rams had a greater contribution to genetic erosion in this population.

The genealogical parameters based on probabilities of gene origin, such as $f_e,$ $f_a,$ $f_{ne}$ have been recognised as important population demographic parameters since their conceptualisation (Boichard et al. 1997) and have been widely used in many studies on genetic diversity in livestock populations (Toro et al. 2014). The parameters $f_e,$ $f_a,$ $f_{ne}$ are very sensitive to the subset of populations used in their estimation, so the ratios $f_e/f,$ $f_e/f_a,$ and $f_{ne}$/$f_e$ are more conclusive. By simultaneously taking into account the estimates of the above parameters, it is obvious that the Istrian sheep has suffered a substantial loss of genetic variability. The $f_e/f$ of 0.37 (rams) and 0.15 (ewes) indicate disequilibrium between founder contributions and evidence the overuse of a small number of animals as parents over many generations. This result was not surprising when considering the size of this closed population and the practices of breeders in their flocks. In many of the flocks encompassed by this study, breeders rarely purchase rams from another flock, thus increasing inbreeding and reducing genetic variability. The $f_e/f$ has been ‘poorly’ reported for sheep populations (e.g. Goyache et al. (2003) reported 0.24 in Xalda, Ghafouri-Kesbi (2010) reported 0.14 in Zandi, and Mokhtari et al. (2014) reported 0.086 in Iran-Black sheep). The most commonly reported measure was $f_e/f_a\mathrm{ }.$ This parameter has been associated with the ‘bottleneck’ effect in wild populations, analogous to the anthropogenic ‘culling’ force in domestic animals. The $f_e/f_a$ approaching to unity represents a desirable state of the population from a conservation point of view and our estimate for this parameter was quite high for both genders (∼1.45). The result obtained indicates a substantial loss of genetic diversity due to a lower than possible number of effectively contributing animals which were used in the reproduction. Significantly lower estimates of $f_e/f_a$ were reported in Morada Nova hair sheep (∼1.0, McManus et al. 2019), Shall (1.001, Hashemi and Hossein-Zadeh 2020), Dorper (1.09, Kiya et al. 2019), Segurena (1.12, Barros et al. 2017), Iran-Black sheep (1.18, Mokhtari et al. 2014), and Moghani breed (1.22, Mokhtari et al. 2015), but there were also reports of considerably higher estimates such as for Baluchi sheep (1.7, Tahmoorespur and Sheikhloo 2011), Xalda (2.02, Goyache et al. 2003), and Kermani sheep (2.07, Mokhtari et al. 2013). Ancestors are considered to have a large marginal contribution to the reference population if their genes were passed through many descendants. According to our estimates, about 50% of the total genetic variation was explained by only 28 rams and 34 ewes. Our pedigree quality control gave substantial credibility to the obtained results, but they undoubtedly deviate from the true parameters due to the lack of pedigree information for some animals. It is ungrateful to compare these results between populations, but it seems that the Istrian sheep lies somewhere between the above-mentioned sheep populations. The $f_{ne}$ which accounts for all random losses of genes during segregations was estimated to be 33 rams and 53 ewes, additionally proving that some sires and dams were used more intensively in reproduction. A similar $f_{ne}$ was reported for Kermani (48, Mokhtari et al. 2013), and higher $f_{ne}$ for Moghani (67, Mokhtari et al. 2015), and Zandi sheep (74, Ghafouri-Kesbi 2010). Simultaneous comparisons (rams and ewes) of estimated $f$ (923, 282), $f_e$ (105, 137), $f_a$ (72, 95), and $f_{ne}$ (33, 53) and their ratios imply substantial loss of founder alleles in this population. However, since founders are always arbitrary, the obtained results must be taken with precaution and should serve as a good approximation (Mrode 2014).

Pedigrees of dairy sheep populations, particularly those bred in low-input systems, are often not as deep as those of dairy cattle which often leads to underestimation of the coefficient of inbreeding ($F$) and inbreeding rate. The $F$ in this population was pretty high at 7.83% and well above the values commonly reported for sheep populations. A high proportion of animals (43%) had an $F$ above 6.25%, exceeding the threshold for ‘passive inbreeding’ advocated by Miglior (2000). The estimated $F$ in the different sheep production types were mostly between 1 and 4 as follows: Zandi and Afshari sheep (1.46% and 1.2%, Ghafouri-Kesbi 2010, 2012), Xalda sheep (1.5%, Goyache et al. 2003), Shall sheep (1.62%, Hashemi and Hossein-Zadeh 2020), Mehraban sheep (1.69%, Yavarifard et al. 2014), Santa Inês sheep (2.33%, Pedrosa et al. 2010), Moghani sheep (2.93%, Hossein-Zadeh 2012), French dairy sheep breeds (2.34–3.11%, Rodriguez-Ramilo et al. 2019), Elsenburg Merino high and low lines (3.9% and 1.6%, Jorgensen et al. 2020), Bharat Merino sheep (2.59%, Gowane et al. 2013), and for seven low numbered French sheep populations (1.9–2.4% Danchin-Burge et al. 2010). Estimates below 1% were reported for Guilan (0.15%, Eteqadi et al. 2014), Dorper (0.32%, Kiya et al. 2019), and Segurena sheep (0.60%, Barros et al. 2017). However, many of the above studies did not report pedigree completeness and criteria for defining the reference population, so inbreeding may have been underestimated due to missing pedigree information.

The Istrian sheep is a small closed population with sufficient ancestral information to detect inbreeding, and the $F$ determined was somewhat to be expected, especially when selection is considered. In small populations, there are generally fewer candidate animals to select from, and breeding without a specifically designed mating plan tends to increase $F,$ especially in trait-selection breeding programs with the mating of superior animals that share more IBD alleles (inherited from common ancestors). The average rate of inbreeding was estimated to be 0.31% which could be threatening to the genetic diversity of the population. It is impossible to completely avoid inbreeding in this low-numbered sheep population, but an exchange of rams between flocks (Goyache et al. 2003) or using a rotational mating scheme (Mokhtari et al. 2013) should be helpful to maintain it at an acceptable level. The above schemes should be helpful not only in maintaining genetic diversity but also in providing genetic links between flocks which is an important prerequisite for a successful between flock BLUP genetic evaluation system (Simm et al. 2001; Eikje et al. 2008).

The average coefficient of co-ancestry can be used as a good proxy for genetic heterogeneity, but it is not as informative as $F$ (e.g. a large number of selected related animals such as sibs produced by ‘re-matings’ do not necessarily imply an increase in homogeneity caused by inbreeding). This is particularly true for optimisation selection programs where the coefficient of co-ancestry plays an important role in generating mating schemes ( Eteqadi et al. 2014; Ghavi Hossein-Zadeh 2016). The estimate of average relatedness (0.0116) was pretty high compared to some other breeds such as Shall (0.000072, Hashemi and Hossein-Zadeh 2020), Guilan (0.000144572, Eteqadi et al. 2014), and Mehraban sheep (∼0.0002, Yavarifard et al. 2014). The reasons are the same as discussed above for $F.$ Selection for dairy yield traits in this population is likely to continue so the further loss of genetic variability is expected. However, the rate of genetic erosion should be slowed, preferably by applying the principles of Optimal Contribution Selection (OCS, Colleau et al. 2004; Hinrichs et al. 2006). The genetic contributions of selected candidates to the next generation should be optimised to balance the benefits and risks of selection by maximising genetic gain and restricting the matting of closely related animals. The earlier the OCS is implemented, the greater the benefit can be expected. Initial steps towards genomic OCS in this population have already been conducted, and the willingness of breeders to participate appears promising for the breed from both perspectives, conservation and selection.

The trait-related genetic parameters play an important role in predicting genetic gain and their estimation is very sensitive to the population genetic background and the statistical model used to estimate them. Accordingly, it is very ungrateful to use the genetic parameters estimated in one population for genetic evaluation (or prediction of genetic gain) in another population, which is why we devoted the second part of the study to this issue. Heritabilities for DMY (0.09), FC (0.03), and PC (0.06) were generally low and to some extent comparable only to estimates reported by Komprej et al. (2009) for three Slovenian dairy breeds, Oravcová et al. (2005) for Improved Valachian, and Othmane et al. (2002) for the Churra breed. In the majority of other studies, heritabilities were considerably higher, and regularly exceeded 0.15 for DMY (e.g. 0.18, El-Saied et al. 1998), 0.20 for FC (e.g. 0.24, Barillet and Boichard 1994), and 0.20 for PC (e.g. 0.20, Hamann et al. 2004). The low heritability in this population likely reflects directional selection that reduces additive genetic variance in the population (so called Bulmer effect, Bulmer 1971). The estimates probably reflect the true parameters in this population, but it is also possible that they were misestimated as a consequence of low connectedness between flocks due to the limited exchange of breeding animals between them. Additive genetic correlations between DMY and FC, DMY and PC, and FC and PC were −0.22, −0.27, and 0.70, respectively. Their direction was generally in agreement with a majority of previous reports, while the magnitude was most consistent with the reports of Brežnik et al. (1999) for Slovenian dairy breeds (–0.29, −0.31, and 0.63) and Oravcová et al. (2005) for Improved Valachian (–0.29, −0.30, and 0.57) and Tsigai breeds (–0.23, −0.27, and 0.58). Although higher milk yield leads to higher protein and fat yield per ewe and still results in higher profitability for breeders, the obtained correlations suggest that further selection only for DMY will definitely ‘dilute’ milk, which might make it unsuitable for cheese manufacturers in Istria. Therefore, we hereby encourage multi-trait approach selection to achieve a balance between milk yield and its’ quality loss. This approach was proposed by the national selection program for this breed but is still rarely applied among the breeders.

The third part of the study was devoted to the impact of inbreeding on dairy traits. Trait-related effects of inbreeding were usually reported for lower fertility, prolificacy, disease resistance, vigour, and survivability in many livestock populations. However, information about the impact of inbreeding on dairy traits in lactating animals has been pretty scarce, particularly for sheep. The results of some previous studies on sheep have led us to hypothesise that inbreeding affects dairy yield traits as well (e.g. the negative correlation between inbreeding and retained heterosis for dairy traits (Murphy et al. 2017) and the negative impact of ewe’s level of inbreeding on the early performance of her lambs (Ercanbrack and Knight 1991; Norberg and Sorensen 2007). By estimating the impact of inbreeding using an animal genetic model, these results represent some of the first direct insights into this issue in dairy sheep populations. Fitting the covariate $F$ in the multivariate statistical model allowed estimation of this effect by simultaneously accounting for all available information affecting phenotype. The statistical significance of the effect examined suggested the potential presence of inbreeding depression, but the very limited magnitude of the effect made the results inconclusive to great extent. These small magnitudes of the effect could be due to a skewed distribution of the $F,$ i.e. the prevalence of animals with low $F.$ The majority of the population (∼ 63%) had $F$ below 12.5%, which could ‘dilute’ the true magnitude of inbreeding, as 10% (‘rule of the thumb’) is considered a benchmark for detecting deleterious effects of inbreeding on reproduction or immune response to disease. For more conclusive results, i.e. to capture the true effect, it would be necessary to have more inbreds with high $F.$ However, we hope that such a scenario will not occur in this population in the future, as we are taking all possible actions to keep $F$ low by following the principles of Optimal Contribution Selection to prevent further genetic erosion and the possible deterioration of ‘fitness’ traits. Further studies on this issue are needed, preferably in already highly inbred populations using a similar statistical approach (animal genetic model) and not by simply regressing $F$ on phenotypes (the most common approach). All parameters estimated in this study suggest that the basic principles of Optimal Contribution Selection should be followed in order to optimally exploit the genetic pool in terms of selection and conservation of this breed.

Conclusion

The analysis of genetic variability and trait-related parameters provided valuable information on past breeding activities in the population of the Istrian sheep. The determined parameters of genetic diversity impose a need to prevent further loss of genetic variability while trying to keep the selection going on. To achieve genetic gain, selection must be based on the estimated breeding values. This is especially important in the presence of the low heritabilities as determined here. Results pertaining to the effect of inbreeding on dairy traits are not entirely conclusive and require additional study prior to any generalisation. In the future, a carefully designed breeding plan needs to be applied to maximise the available genetic pool of this breed, preferably following the basic principles of Optimal Contribution Selection.

Acknowledgments

The authors thank the Croatian Science Foundation for its financial support. The authors also acknowledge the support from Istrian sheep breeders, the Croatian Agency for Agriculture and Food, and the Ministry of Agriculture.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All data analysed during this study are included in this published article. The data reported in this study are available from the corresponding author, Marija Špehar, PhD.

Additional information

Funding

The study was funded by the Croatian Science Foundation (Genomic characterisation, preservation and optimum contribution selection of Croatian dairy sheep, OPTI-SHEEP), grant number IP-2019-04-3559.