Skip to main content
  • Research Article
  • Open access
  • Published:

A comprehensive genome-wide scan detects genomic regions related to local adaptation and climate resilience in Mediterranean domestic sheep



The management of farm animal genetic resources and the adaptation of animals to climate change will probably have major effects on the long-term sustainability of the livestock sector. Genomic data harbour useful relevant information that needs to be harnessed for effectively managing genetic resources. In this paper, we report the genome characterization of the highly productive Mediterranean Chios dairy sheep and focus on genetic diversity measures related with local adaptation and selection and the genetic architecture of animal resilience to weather fluctuations as a novel adaptative trait linked to climate change.


We detected runs of homozygosity (ROH) and heterozygosity (ROHet) that revealed multiple highly homozygous and heterozygous hotspots across the Chios sheep genome. A particularly highly homozygous region was identified on chromosome 13 as a candidate of directional genetic selection associated with milk traits, which includes annotated genes that were previously shown to be linked to local adaptation to harsh environmental conditions. Favourable heterozygosity related with a potentially protective role against livestock diseases and enhanced overall fitness was revealed in heterozygous-rich regions on sheep chromosomes 3, 10, 13 and 19. Furthermore, genomic analyses were conducted on sheep resilience phenotypes that display changes in milk production in response to weather variation. Sheep resilience to heat stress was a significantly heritable trait (h2 = 0.26) and genetically antagonistic to milk production. Genome-wide association and regional heritability mapping analyses revealed novel genomic markers and regions on chromosome 5 that were significantly associated with sheep resilience to climate change. Subsequently, an annotation analysis detected a set of genes on chromosome 5 that were associated with olfactory receptor complexes that could participate in heat stress mitigation through changes in respiration rate and respiratory evaporation. Other genes were grouped in previously reported biological processes relevant to livestock heat dissipation, including stress and immune response.


Our results may contribute to the optimal management of sheep genetic resources and inform modern selective breeding programmes that aim at mitigating future environmental challenges towards sustainable farming, while better balancing animal adaptation and productivity. Our results are directly relevant to the studied breed and the respective environmental conditions; however, the methodology may be extended to other livestock species of interest.


Intensive genetic selection of farm animals combined with environmental changes raises concerns about the preservation of livestock genetic diversity. Modern programmes of selective breeding must integrate the effective management of genetic resources and the adaptation of animals to present and future challenges [1]. The integration of methods to mitigate genetic erosion and comprehensively evaluate animal performance through the characterization of adaptive traits that are relevant to climate change adaptation [2] would facilitate the sustainable management of livestock biodiversity [3].

Assessment of genome diversity, by estimating homozygosity and heterozygosity levels and genomic inbreeding, is of paramount importance for the management and conservation of genetic resources [4]. The availability of medium- and high-density genome-wide DNA arrays for commercial livestock species [5, 6] enables the estimation of multiple genomic diversity indicators. In particular, the effective quantification of individual homozygosity levels based on the detection of runs of homozygosity (ROH) [7, 8] could be used to assess and manage diversity and prevent inbreeding depression [9]. Estimation of inbreeding coefficients based on ROH is considered as a powerful approach that is widely used in biodiversity studies [10,11,12], and offers the possibility to distinguish between ancient and recent inbreeding [13]. In addition, the detection of ROH hotspots (also referred to as ROH islands) can identify suggestive genomic selected regions that harbour genes associated with animal domestication [14], breed formation [15], traits of economic importance [16], and local adaptation to new environments [12]. Furthermore, heterozygosity levels are commonly estimated as a genetic diversity indicator. Recently, it has been proposed that highly heterozygous regions throughout the genome, known as heterozygosity-rich regions (ROHet) or ROHet islands [17], are strongly associated with animal fitness and survival, and heterotic balancing selection processes [17,18,19].

The study of the genetic architecture of adaptive traits that are relevant to animal performance in response to changing environmental conditions, such as weather fluctuations, could also contribute to the characterization of animal genetic resources [20], and at the same time, provide essential information for future selective breeding strategies to mitigate the effects of climate change on the respective breeding goals [21, 22]. Genome-wide association studies on heritable adaptive traits have recently focused on animal production resilience to climate change in dairy cattle [23] and dairy goats [24], and have revealed potential candidate genes that are involved in animal stress and metabolic responses. Simulation studies have demonstrated that the integration of resilience traits in selective breeding programmes could enhance long-term sustainability of animal performance [22].

The present study focuses on the genomic characterization of the Chios dairy sheep, which is one of the most productive sheep breeds worldwide. Our specific objectives were to: (i) derive and assess multiple genetic diversity indicators including levels of homozygosity and heterozygosity, and genomic inbreeding estimates, (ii) detect genomic regions associated with local adaptation and selection based on ROH and ROHet hotspots, and (iii) characterize an adaptive trait relevant to climate change defined as the performance resilience of animals to weather volatility.


DNA extraction and genotyping

Total genomic DNA was extracted from blood samples of 600 Chios female sheep (milking ewes) from three farms in Northern Greece, according to the standard protocol of the QIAamp DNA Mini and Blood Mini kit (QIAGEN, USA), following the manufacturer’s recommendations with minor modifications. DNA samples were quantified using a Nanodrop spectrophotometer and stored at − 20 °C until use. Subsequently, DNA samples were oven-dried and 500 ng per sample were transferred to 96-well plates. All DNA samples were genotyped with the OvineSNP50 Genotyping BeadChip that contains 54,124 single nucleotide polymorphisms (SNPs) (Illumina, Inc., U.S.) and was developed in collaboration with the International Sheep Genomics Consortium [25] which included 23 Chios sheep.

Quality control of genotypes was performed by setting sample and marker call rate thresholds at 90% using the PLINK 1.07 software [26]. In addition, all SNPs on the sex chromosomes were removed resulting in a final dataset of 538 individuals and 51,124 SNPs spread across the 26 ovine autosomes. SNP positions were assigned according to the Oar_v4.0 sheep genome assembly [27].

Genome-wide characterization

First, genotypes were analysed to derive observed (Ho) and expected (He) heterozygosity levels, and the inbreeding coefficient (F) based on the difference between the observed and expected number of homozygous genotypes under Hardy–Weinberg equilibrium, which were calculated from the allele frequencies in the sample and the total number of markers. These analyses were conducted with the PLINK 1.07 software [26].

Subsequently, ROH and ROHet were detected by implementing the R package detectRUNS [28] and applying the sliding-window approach for ROH and the consecutive method for ROHet detection [29]. The following input parameters were set to detect ROH based on a previous study for ROH identification that used medium-density SNP data [30]: (i) a sliding window size of 50 SNPs across the genome; (ii) a proportion of homozygous overlapping windows of 0.05; (iii) for each ROH, a minimum number of 20 consecutive SNPs and a minimum length of 250 kb; (iv) a maximum gap of 250 kb between consecutive homozygous SNPs; and (v) a maximum number of one missing call and one heterozygous genotype allowed in a ROH. Subsequently, highly homozygous regions, defined as ROH islands, were identified by extracting the 99.9 quantile of the frequency distribution of SNPs within a ROH.

The genomic inbreeding coefficient based on ROH (\({F}_{ROH}\)) was also estimated for each individual, as described by McQuillan et al. [31]:

$${F}_{ROH}=\frac{\sum {L}_{ROH}}{{L}_{aut}},$$

where \(\sum {L}_{ROH}\) is the total length of all ROH in each individual genome and \({L}_{aut}\) is the total length of the autosomal genome according to the DNA array used for genotyping.

Contrary to ROH, to date the optimal parameters for ROHet detection have not been extensively studied. In the current study, we set the following input parameters: (i) for each ROHet, a minimum number of 15 consecutive SNPs and a minimum length of 250 kb; (ii) a maximum gap of 250 kb between consecutive homozygous SNPs; and (iii) a maximum number of two missing calls and three homozygous genotypes. Highly heterozygous genomic regions, defined as ROHet islands, were identified based on the frequency distribution of SNPs within a ROHet, as described for ROH islands above.

Characterization of climate resilience

For the subsequent analyses, an additional edit was applied to the genotyping data, where SNPs with a minor allele frequency lower that 0.01 were removed. Thus, the remaining analyses were based on the same 538 individuals and the remaining 47,605 SNPs across the 26 ovine autosomes.

Genotype principal component analysis

A principal component analysis was performed on the animal genotypes to explore the potential population structure that should be accounted for in the following downstream genomic analyses to avoid possible inflation effects. The GEMMA software [32] was used for this purpose and the decomposed relatedness matrix between individuals was visualized.

Animal phenotypes

Milk yield records were available for all genotyped animals during the 2003–2018 period, during which animals had multiple lambings (births of lambs) followed by lactation periods (milking periods) lasting on average five months. Milk yield was recorded on each individual on a specific date once per month within a lactation. Thus, multiple milk yield records were available for each genotyped animal. This dataset comprised 10,348 records on 538 individuals.

Corresponding weather data were obtained from the Hellenic National Meteorological Service pertaining to the weather stations that were closest to the three farms on which animals were sampled. Average daily air temperature data were extracted for the dates corresponding to the available milk yield records.

Individual resilience phenotypes, reflecting a change in milk production in response to temperature variability, were derived based on a methodology that is described in detail in [33]. Briefly, reaction norm functions were fitted to milk yield records using a random regression model that included the corresponding air temperature measurement and accounted for the fixed effects of farm, lactation number, calendar year and month of lambing, and the number of days from lambing when milk yield was recorded. Resilience phenotypes were essentially the resulting slopes of individual reaction norm curves at different temperature levels. For the purposes of the present study, we focused on temperatures of 10 °C and 25 °C, which represent thresholds of cold and heat stress, respectively, for sheep raised in this geographic region. Animal resilience reflected the changes in milk yield in response to changes in temperature under cold and heat stress conditions. The above analyses were conducted using the BLUPF90 software [34] in RStudio with R.3.6.1 [35].

The newly derived resilience phenotypes were studied together with animal lifetime milk production, estimated from daily milk yield records according to previously described formula [36], and length of productive life, defined as the total number of days a ewe was milked.

Estimation of genomic parameters and breeding values

Variance components of animal resilience to temperature variability in cold and hot weather conditions were estimated to implement the following mixed linear model:


where \(\mathbf{y}\) is the vector of animal phenotypes, \({\varvec{\upalpha}}\) is the vector of associated fixed effects, \(\mathbf{u}\) is the vector of random polygenic (additive genetic) effects, which is distributed as a multivariate normal distribution \(MVN({\bf{0}},{V}_{g}\mathbf{G})\) with \(\mathbf{G}\) being the genomic relationship matrix and \({V}_{g}\) the genomic variance of the trait, \(\mathbf{W}\) and \(\mathbf{Z}\) are the corresponding design matrices, and \({\varvec{\upepsilon}}\) is the vector of random residual effect. Fixed effects included in Model (1) were: the farm on which the animals were raised, calendar year and month of first lambing, total number of lactations, and length of productive life; the first three principal components extracted from the population structure analysis described above were also included as covariates.

In separate analyses, variance components were also estimated for lifetime milk production and length of productive life using Model (1), with the addition of the latter as a covariate in the analysis of the former, and vice versa. Thus, results pertained to lifetime milk production adjusted for length of productive life, which reveals the true producing capacity of the animals. Furthermore, length of productive life was adjusted for the amount of milk produced, which then represents a functional longevity proxy. This trait reflects the capacity of an animal to remain on the farm independently of its production level, implying good health and fitness.

Genomic relationships between animals were estimated based on their SNP genotypes using the Regional Heritability Advanced Complex Trait Analysis (REACTA) software [37] and were incorporated in the analyses above. In each case, trait heritability was derived as the ratio of the estimated additive genetic variance to the total phenotypic variance.

Phenotypic and genomic correlations between the studied traits were subsequently estimated in a series of bivariate analyses with the above-described models. Statistical significance of all estimates was assessed using the two-tailed Student’s t-distribution.

The same model was also used to derive the genomic estimated breeding values (GEBV) of each individual for the studied traits. For resilience traits, these values represent the inherent propensity of animals to maintain milk production levels regardless of changing weather conditions. The respective GEBV accuracy was estimated as:

$$A\mathrm{c}curacy=\sqrt{1-\frac{PEV}{{\sigma }_{G}^{2}}},$$

where \(PEV\) is the prediction error variance of the GEBV and \({\sigma }_{G}^{2}\) is the trait additive genetic variance. In this context, accuracy of the GEBV reflects the correlation of the GEBV with the true unknown genetic value of an animal. All these analyses were conducted using the ASReml 4.1 software [38].

Genome-wide association study

Animal phenotypes and genotypes were jointly analysed to identify SNPs associated with animal resilience to hot and cold weather, lifetime milk yield, and length of productive life. The same effects as in the above-described Model (1) were fitted, with the addition of the term \(\mathbf{x}{\varvec{\upbeta}}\), with \(\mathbf{x}\) being the vector of SNP genotypes and \({\varvec{\upbeta}}\) their associated effects. Analyses were conducted using the GEMMA software [32]. A Bonferroni correction was applied for multiple testing to determine a genome-wide (P < 0.05) and suggestive (one false positive per genome scan) significance threshold, resulting in final threshold values of P < 1.05E−06 and P < 2.10E−05, respectively, corresponding to −log10 (P) of 5.98 and 4.68, respectively. In order to account for further cryptic population structure that may inflate the test statistics, the inflation factor λ was used for correction [39].

Regional heritability mapping

Regional heritability mapping (RHM) was implemented to identify potential genomic regions associated with animal resilience and the other traits by scanning windows across the whole genome. A sliding-window approach was used with genomic regions of 100 consecutive SNPs (several window sizes were tested to evaluate the optimal region size) and overlaps of 30 SNPs defined along each autosome. The REACTA software [37] was used with the following model:


where \(\mathbf{y}\) is the vector of individual animal phenotypes, \({{\mathbf{u}}}_{(i)}\) is the vector of the random effect of region \(i\), distributed as \(MVN[\bf{0},{\mathrm{V}}_{\mathrm{g}(\it {i}\mathrm)}{\mathbf{G}}_{(\it{i}\mathrm)}]\) (with \({\mathrm{V}}_{\mathrm{g}(i)}\) and \({\mathbf{G}}_{(i)}\) being the genomic variance and the genomic relationship matrix corresponding to the SNPs in region \(i\), respectively) and \({\mathbf{u}}_{(-i)}\) is the random polygenic effect excluding region \(i\), distributed as \(MVN[\bf{0},{\mathrm{V}}_{{\mathrm{g}}({-}\it{i}\mathrm)}{{\mathbf{G}}}_{({-}\it{i}\mathrm)}]\) (with \({\mathrm{V}}_{\mathrm{g}(-i)}\) and \({\mathbf{G}}_{(-i)}\) being the genetic variance and the genomic relationship matrix corresponding to all SNPs other than those in region \(i\), respectively); all other effects are as described for Model (1).

The likelihood ratio test (LRT) was used to assess the significance of the effect of each genomic region. In total, 681 regions were tested across the genome resulting in a genome-wide significance threshold (P < 0.05) after Bonferroni correction for multiple testing corresponding to P < 7.34E−05 with −log10 (P) of 4.13 and a suggestive significance threshold (one false positive per genome scan) defined at P < 1.47E−03 with −log10 (P) of 2.83.

Gene and functional annotation analysis

The extent of linkage disequilibrium (LD) was determined in regions where (i) ROH or ROHet islands were detected, and (ii) genome-wide and suggestive significant SNPs were identified in genomic and RHM analyses. Pairwise LD estimation (r2) and visualization were carried out using the PLINK v1.07 [26] and Haploview 4.2 [40] software. Then, gene annotation was performed using the Variant Effect Predictor software [41] within distances of ± 250 kb and 1 Mb according to pairwise LD results of ROH and ROHet, and genomic and RHM analyses, respectively. Functional annotation was performed for all the identified genes within the respective regions using the UniProt database [42]. Functional enrichment analyses were conducted using the gProfiler genomic tool [43] and the DAVID software [44] to detect candidate annotated genes involved in biological processes with relevance to animal resilience to climate change and the other biological functions.


Genome-wide characterization

The average observed (Ho) and expected (He) heterozygosity estimates were 0.314 ± 0.167 and 0.355 ± 0.201, respectively, revealing quite high levels of polymorphism within Chios sheep breed. Average ROH-based inbreeding (\({F}_{ROH}\)) was 0.1133 and individual animal estimates ranged from 0.0314 to 0.3013. Estimates of the inbreeding coefficient (F) based on the difference between the observed and expected number of homozygous genotypes ranged from − 0.1123 to 0.2465, with an average of 0.005. Inbreeding levels were also estimated for different lengths of ROH segments and were lower for longer ROH (Table 1). Short ROH display ancient inbreeding and long ROH represent more recent inbreeding events [45]. High Pearson correlations were estimated between F and \({F}_{ROH}\) (0.97), and between coefficients of recent inbreeding (0.98), whereas lower estimates were obtained between recent and ancient inbreeding (0.78–0.80) (see Additional file 1: Fig. S1). Figure 1 illustrates the individual ROH-based inbreeding levels of Chios sheep according to farm of origin and shows that the population in farm C is less inbred, possibly implying a more intensive human-mediated directional selection on the other two farms.

Table 1 Average genomic (F) and ROH-based (FROH) inbreeding coefficients
Fig. 1
figure 1

Individual run of homozygosity (ROH)-based inbreeding levels (FROH) in Chios sheep. A, B, C correspond to the farm of origin of the animals

In total, 23,725 ROH were identified throughout the genome. The average number of ROH per animal was 44.1 and the average ROH length detected across all the autosomes was 4.92 Mb. The shortest ROH (736.2 kb) was detected on chromosome 14 and the longest (97.57 Mb) on chromosome 2. The highest individual sum of ROH reached 737.52 Mb with the corresponding maximum \({F}_{ROH}\) of 0.3013. The numbers of ROH identified by size class (length of ROH) are in Table 2.

Table 2 Number of identified ROH and ROHet per size class

Furthermore, 15,218 ROHet were identified in the studied sample of animals. The average ROHet length detected across all the autosomes was 719.92 kb; the shortest was detected on chromosome 1 (376.04 kb) and the longest on chromosome 10 (1748.15 kb). The numbers of ROHet by size class (length of ROHet) are in Table 2.

Figure 2 illustrates the frequency of SNPs within a ROH or a ROHet across the autosomes, revealing homozygosity-rich and heterozygosity-rich genomic regions, respectively. These regions, termed ROH and ROHet islands, were identified by obtaining the 99.9 quantile of the frequency distributions of SNPs within a ROH or ROHet, corresponding to a frequency threshold of 0.61 and 0.32, respectively. Five such genomic regions were identified: (i) one ROH island representing a genomic region of 2.89 Mb on chromosome 13, where 55 SNPs and 57 genes are localized; and four ROHet islands, i.e. (ii) one (681.96 kb) on chromosome 3 including 15 SNPs, (iii) one on chromosome 10 (489.35 kb) including 12 SNPs, (iv) one on chromosome 13 (808.68 kb) including 16 SNPs and 25 annotated genes, and (v) one on chromosome 19 (479.36 kb) including nine SNPs and seven closely located genes. High LD levels were observed in all these chromosomal regions (see Additional file 2: Table S1).

Fig. 2
figure 2

Manhattan plots showing the frequency of a SNP within homozygous and heterozygous segments in Chios sheep. a Runs of homozygosity (ROH), and b runs of heterozygosity (ROHet). Grey lines indicate the 99.9 quantile thresholds, defining the respective ROH/ROHet islands

Functional annotation analyses—genomic characterization

Based on the average pairwise LD (r2) and average distance between the SNPs detected in the ROH and ROHet islands, a SNP annotation analysis was performed in 200–250 kb-regions upstream and downstream of the SNPs. Ninety-nine genes were located in the annotated regions of the respective ROH and ROHet islands (see Additional file 2: Table S1). Three of the candidate annotated genes located in the ROH island on chromosome 13 (BMP7, ACSS2, and EPB41L1) are known to be associated with animal traits, such as litter size in sheep [46] and milk quality in sheep and cattle [32, 47, 48], while the other genes (CTCFL, PCK1, RAE1, and BMP7) are reported to be involved in local adaptation processes [49]. Moreover, the annotated genes located in highly heterozygous regions (ROHet islands) on chromosomes 13 and 19, respectively, are potentially linked with animal immune response (ROMO1), protection against ruminant disorders (LOC101117181, and SUMF1), and cheese making traits (ITPR1) [50,51,52,53]. In addition, genes related to sheep domestication (ITPR1 and SUMF1) [54] were found in the ROHet island on chromosome 19, implying positive selection on heterozygous genotypes. All the candidate genes that were identified in ROH and ROHet islands and that are associated with livestock performance and disease and with adaptation processes are summarized in Additional file 2: Table S1. All annotated genes and their molecular and functional characterization are also summarized in Additional file 3: Table S2.

Characterization of climate resilience

Population structure

The results of the principal component analysis that was carried out on the genomic relatedness matrix between individuals are illustrated in Fig. 3. The first three principal components accounted for almost 11% of the total variation (Fig. 3a). A strong population structure was revealed by the first two principal components and attributed to the farm of origin of animals (Fig. 3b). The latter may be explained by targeted genetic selection and mating practices on the different farms. A similar structure was revealed by the first and third principal components. Therefore, the first three principal components were fitted as covariates, together with the genomic relatedness matrix, in the following downstream genomic analyses to account for population structure.

Fig. 3
figure 3

Principal component analysis (PCA) of Chios sheep. a Proportion of variation (%) corresponding to each principal component of the decomposed genomic relationship matrix; and b First and second principal components showing the population stratification attributed to the farm origin of individual sheep

Animal resilience phenotypes

Table 3 summarizes the statistics for the traits under study in the Chios dairy sheep dataset. The mean (SD) of average daily temperature was 15.51 °C (8.45). Figure 4 illustrates the average Chios sheep performance resilience to temperature fluctuations, which reflects the change in population daily milk yield in response to temperature volatility. Examples of individual animal performance resilience phenotypes represented by the slopes of reaction norm curves (see Additional file 4: Figure S2) reflect deviations from the population curve shown in Fig. 4.

Table 3 Descriptive statistics (mean and standard deviation) of the traits under study in the Chios dairy sheep dataset
Fig. 4
figure 4

Chios sheep milk performance resilience to air temperature variation. The curve represents the average change in daily milk yield (DMY) in response to temperature variation (Tavg, °C)

Indicative resilience phenotypes to cold (10 °C) and hot (25 °C) conditions were selected for further study. Positive or negative values of individual slopes represent milk yield changes in response to ambient temperature change, whereas zero (flat) slopes represent individuals whose milk production remains unaffected by temperature fluctuations.

Genomic parameters of the traits under study

Trait heritability and correlation estimates are in Table 4. Heritability estimates of animal resilience to hot weather and lifetime milk yield were moderate and significantly higher than 0 (h2 = 0.20 and 0.26, SE = 0.095 and 0.009, respectively). Resilience to cold temperature and length of productive life were lowly heritable with estimates that did not reach statistical significance. The two animal resilience traits were strongly positively correlated at both the phenotypic and genetic levels. Phenotypic and genetic correlations of animal resilience under hot conditions with lifetime milk yield were negative, implying an antagonistic relationship between sheep production and fitness, while those between resilience to hot weather and length of productive life were positive as they both represent fitness traits. All other correlations were not statistically different from 0.

Table 4 Heritability (diagonal) and correlation (above: genomic, below: phenotypic) estimates between the traits under study with standard errors in parentheses

Genomic breeding values

The comparison between individual GEBV estimated for the different studied traits is illustrated in Fig. 5. Notably, positive breeding values for resilience to cold and hot weather correspond to low genetic merit for lifetime milk yield and vice versa, implying a likely genetic antagonism for these traits, which is consistent with the negative genetic correlations in Table 4. The opposite is true for resilience to hot weather and length of productive life. In addition, high breeding values for length of productive life correspond to lower genetic merit for lifetime milk yield in accordance with the estimated negative correlation between these traits (Table 4).

Fig. 5
figure 5

Comparison of estimated genomic breeding values (GEBV) of Chios sheep resilience to cold (Tavg10) and hot (Tavg25) weather conditions, lifetime milk production (tmilk), and length of productive life (tdim). a Resilience to cold and b hot weather compared to lifetime milk production respectively, c comparison between resilience to cold and hot weather conditions, d comparison of resilience to cold, e hot weather with length of productive life respectively, and f comparison between lifetime milk production and length of productive life

Average GEBV accuracy estimates were 0.47, 0.20, 0.56 and 0.34 for animal resilience to hot and cold weather, lifetime milk yield, and length of productive life, respectively. The magnitude of these parameters reflects the corresponding trait heritability estimates summarized in Table 4. Higher heritability estimates are associated with more accurate GEBV.

Genome-wide association study

Results from the GWAS of animal resilience under hot and cold conditions, lifetime milk yield, and length of productive life are illustrated in Fig. 6. The corresponding Quantile–Quantile (Q–Q) plots are provided in Additional file 5: Figure S3 and show no significant inflation. The results reveal a largely polygenic mode of inheritance for these traits, with specific genomic regions of interest located on chromosomes 5 and 19 for resilience to hot weather and lifetime milk performance.

Fig. 6
figure 6

Genome-wide association results for Chios sheep resilience to temperature change and performance-related traits. a Resilience to hot weather conditions, b resilience to cold weather conditions, c lifetime milk yield, and d length of productive life. Red and blue lines indicate the genome-wide and suggestive (one false positive per genome scan) significance thresholds, respectively

One SNP approached the genome-wide significance threshold (P = 1.15E−06) and five genome-wide suggestive significant SNPs were also detected for resilience to hot weather. Interestingly, one of the latter on chromosome 19 also had a genome-wide significant association with lifetime milk yield with a substitution effect of opposite sign compared to resilience (Table 5), in concordance with the negative genetic correlation between the two traits (Table 4; Fig. 5). No significant SNP associations were detected for resilience to cold weather and length of productive life, consistently with the non-significant heritability estimated for these traits (Table 4).

Table 5 Genome-widea and suggestive significant SNPs associated with animal resilience to temperature change in hot weather and lifetime milk yield

Regional heritability mapping

Two neighbouring genomic regions were identified on chromosome 5 that exceeded the genome-wide and suggestive significance threshold (Fig. 7 and Table 6) for resilience to temperature change under hot weather conditions. Four of the five significant SNPs detected by GWAS were located within one region and the fifth was located in the other. Respective regional heritability estimates ranged from 0.06 to 0.07 (Table 6). Furthermore, a region exceeding the suggestive threshold was detected on chromosome 19, including SNP s14444.1, which was suggestive significant in the GWAS.

Fig. 7
figure 7

Regional heritability mapping results for Chios sheep resilience to temperature change and performance-related traits. a Resilience to hot weather conditions, b resilience to cold weather conditions, c lifetime milk yield, and d length of productive life. Red and blue lines indicate the genome-wide and suggestive (one false positive per genome scan) significance thresholds, respectively

Table 6 Regional heritability (h2reg) for the detected genomic regions significantly associated with animal resilience to temperature change in hot weather, lifetime milk yield, and length of productive life

A genomic region located on chromosome 13 was detected for lifetime milk yield exceeding the genome-wide significance threshold, although this region did not include any SNPs identified as significant in the GWAS. A suggestive significant genomic region was detected on chromosome 1 for length of productive life. RHM analysis of animal resilience under cold temperatures did not reveal any genome-wide or suggestive significant associations.

Functional annotation analyses—climate resilience

The LD structure of the region on chromosome 5 where the significant SNPs were identified in the genomic analysis of animal resilience to temperature change in hot weather was estimated (see Additional file 6: Figure S4). According to the results of pairwise LD, a SNP annotation analysis was performed within 1-Mb regions upstream and downstream of each SNP, and 94 candidate genes were detected in this region (see Additional file 7: Table S3). Then, this gene list was used as input to the functional enrichment analysis. The resulting protein interaction network includes all the identified candidate genes (see Additional file 8: Figure S5). Some of the identified biological processes are relevant to adaptation to ambient temperature fluctuations and to olfactory processes (see Additional file 9: Figure S6). Indicatively, the olfactory receptor 56-like and Olfr56 genes are located 97,975 and 167.41 kb upstream of the suggestive significant SNP OAR5_41355579.1, whereas the olfactory receptor 7A17-like gene is located 212.78 kb downstream. In addition, several other annotation clusters were detected (see Additional file 10: Table S4) including genes that participate in biological processes that are related with stress response, such as the glucose, energy and protein metabolism (B4GALT7, HK3, FAF2, MGAT1, and LMAN2), and immune response (BTNL9, TRIM41, TRIM52, FAF2, and PRR7) [55, 56].


The present study was set out to provide a first insight into the genome of the highly productive Chios dairy sheep breed, including an assessment of its genetic diversity, identification of regions related with adaptation and selection, and dissection of the genomic architecture of an adaptive trait linked to animal resilience to climate change.

We derived multiple genome-wide indicators of the Chios genome diversity. Average observed and expected heterozygosity levels (Ho = 0.31, He = 0.35) were consistent with previously reported estimates for domestic sheep breeds worldwide [57] and particularly with other dairy sheep breeds raised in the Mediterranean region (He = 0.38–0.39) [12, 58]. By extension, this applies to expected and observed homozygosity levels which were used to derive the genomic inbreeding coefficient (F) shown in Table 1. Moreover, estimated average ROH-based inbreeding levels (\({F}_{ROH}\)) were moderate (Table 1 and Fig. 1) and within the range of previous ROH-based estimates in domestic dairy sheep (0.04–0.15) [12, 59, 60]. It should be noted here that these inbreeding coefficients collectively reflect genetic diversity levels in the population, but individually have different definitions and scales of measurement. Thus, FROH is an identity by descent probability ranging from 0 to 1, whereas, F is estimated based on allele frequencies resulting in both positive and negative values; this may explain the near zero mean value of the latter in our data (Table 1).

Inbreeding levels in farmed livestock are often attributed to selective pressure towards enhanced productivity resulting in the likely selection of and mating between genetically-related individuals. Long-term unmonitored selection towards the same direction, for example increased productivity, may result in increased homozygosity and compromise animal fitness [61]. In the present study, the strong population structure attributed to farm origin (Fig. 3b) implies mainly within-farm animal genetic selection and reproduction.

Patterns of homozygosity were observed according to different clusters of ROH by size class (Table 2 and Fig. 2a). Short to medium size ROH (0–6 Mb) constitute the majority of the identified homozygous segments (78%) and are indicative of ancient inbreeding events [31] or random recombination events in the respective genomic regions [62]. Long ROH were also present (Table 2), which indicate recent common ancestry due to a low probability of recombination events and selective pressure [14, 62].

We detected a genomic region with increased homozygosity, defined as a ROH island, spanning 2.89 Mb on chromosome 13 of the Chios sheep genome (Fig. 2a) (see Additional file 2: Table S1). This ROH island coincides with a genomic region previously identified as a signature of selection in European dairy sheep breeds [63], where quantitative trait loci (QTL) associated with milk quality traits are localized. These results suggest that a common selective pressure occurred independently in different domestic dairy sheep breeds.

Contrary to ROH islands, heterozygosity-rich regions (ROHet islands), which may be associated with balancing selection favouring animal fitness, are poorly characterized in livestock and, to our knowledge, have not been previously examined in Mediterranean sheep [64]. In this study, we observed patterns of clustered heterozygosity that reveal highly polymorphic regions of different lengths across the genome (Table 2 and Fig. 2b). Regions of 500 to 1000 kb constituted the majority of such heterozygosity-rich segments (91%). The maximum length of an identified ROHet island was 2 Mb, which is consistent with previously reported results for cattle [10] and horse [19] breeds. It is worth mentioning that the detection of ROHet and the respective selected parameters have not been extensively studied yet. Indeed, to date only one study including a sensitivity analysis on specific parameters, such as the numbers of missing and homozygous genotypes allowed within a ROHet segment, has been performed on cattle, suggesting that these regions could be characterized as heterozygosity-rich regions with expected homozygous loci distributed within the segments [10].

We detected four ROHet islands on ovine chromosomes 3, 10, 13 and 19, providing signals of favoured heterozygosity related to protection against livestock diseases, immune response, and local adaptation. Particularly, the heterozygous-rich segment on chromosome 10 (489.35 kb) is mapped to a selection sweep region that was previously identified in Egyptian sheep and includes a region of conserved synteny between Egyptian sheep and goat breeds adapted to hot arid environments [65]. Adaptation to such extreme local environments may involve the presence of multiple heterozygous genotypes associated with improved animal fitness in the respective conditions. Furthermore, the ROHet island located on chromosome 13 includes a region that was previously characterized as a signature of selection in local Russian sheep [66] and is related to mammary gland function and milk fatty acid composition [48]. The observed heterozygosity in this region may reflect a favourable impact on udder health compared to homozygous recessive genotypes.

Functional annotation analysis of the detected ROH island on chromosome 13 revealed a set of genes previously characterized as signatures of selection for local adaptation to harsh environmental conditions (CTCFL, PCK1, RAE1, BMP7, and SPO11) in sheep [49], and animal traits of economic interest related with production and mammary gland development and differentiation (NFS1 and TFAP2C) [63, 67]. Interestingly, the BMP7 and TFAP2C genes are also associated with lifetime milk production in goats [24]. This highly homozygous chromosomal segment may indicate a region of positive selection possibly attributed to selective breeding for increased milk production. Future comparative ROH studies on other dairy sheep breeds raised in similar environments could clarify the importance of this genomic region. Furthermore, annotated genes located within or close to the identified ROHet islands were mostly associated with immune response, protection against livestock disorders, and other animal performance traits (see Additional file 2: Table S1). The detected ROHet islands imply that heterozygosity is favoured for these genomic regions involved in common livestock diseases. For example, QTL related to susceptibility to ovine transmissible spongiform encephalopathy [68] and bovine spongiform encephalopathy [69] are located close to a ROHet island detected in the present study. This ROHet island also harbours genes that are associated with other disorders (LOC101117181 and CEP250), metabolic pathways (EIF6), and milk-related traits (ACSS2, and EPB41L1), which suggest a multiple protective role. Highly heterozygous regions may also be associated with natural selection contributing to overall fitness in contrast with the ROH islands where directional genetic selection possibly leads to high levels of homozygosity.

To further investigate Chios sheep adaptiveness to changing environmental conditions, we evaluated the performance resilience of animals to weather fluctuations (Fig. 4) with a view to possibly include adaptive traits in future selective breeding schemes. For this purpose, we derived slopes of reaction norm functions fitted to random regression models, which reflect changes in individual milk production in response to changing daily temperatures (see Additional file 4: Fig. S2). The effect of heat stress is of particular concern to Mediterranean sheep production [56, 70, 71]. Notably, current and projected scenarios of climate change in this geographic area are suggestive of rising air temperature and increased weather volatility [72]. The heritability estimate of Chios sheep resilience to temperature change in heat stress conditions was significantly higher than zero (Table 4). Although these genomic estimates were based on a relatively small number of genotyped individuals, they are consistent with previous population estimates derived from a pedigree analysis of approximately 37,000 animals [33]. These results confirm that this trait could be improved through genetic selection. However, caution needs to be exercised when combining resilience with other animal traits in selective breeding. Phenotypic and genetic correlations between traits will determine the optimal way of combining them towards holistic animal improvement. In the present study, an antagonistic genetic correlation of resilience to hot weather with lifetime milk production was estimated (Table 4), which is consistent with previous estimates between animal resilience and milk production using pedigree data [33, 73]. This genetic antagonism was also evident in the estimated genomic breeding values for the two traits (Fig. 5), with sheep that are genetically predisposed for the highest milk production having negative values for resilience and length of productive life. This means that if selective breeding continues to focus on improving milk production, resilience to hot weather will decrease, as exemplified by reduced milk yield under heat stress conditions. Length of productive life that reflects animal longevity, would also decrease. Nevertheless, the results on length of productive life should be viewed with caution because of the lack of statistically detectable genetic variation for this trait in the studied population. These results imply the need for a multi-trait breeding index to underpin improvement in both animal production and fitness traits, where the latter are manifested in milk yield being unaffected by weather volatility and long lifespan. The feasibility of developing such indices in selective breeding given the appropriate input parameters has been documented [22].

Our genome-wide association (GWAS) and regional heritability mapping (RHM) analyses identified several SNPs, regions, genes, and molecular pathways associated with Chios sheep resilience under hot weather conditions (Figs. 6a and 7a, and Tables 5 and 6). These results may be used to enhance the accuracy of derived genomic breeding values in a selective breeding programme. In addition, one common SNP was associated with both resilience to hot weather and lifetime milk production and had an opposite substitution effect (Table 5), which is consistent with the above-mentioned genetic antagonism between the two traits. This information should be considered when combining these two traits in a breeding index for selection purposes. It is worth mentioning that very few genomic studies on animal resilience to climate change have been conducted [23, 24] and none of these reported regional heritability mapping results. Regional heritability mapping has been proposed as a complementary approach to detect genomic regions harbouring causative alleles of small effect that contribute to trait variation and could not be captured individually by GWAS [74]. Characteristically, our RHM results detected genomic regions on chromosomes 5 and 19 that are related to animal resilience to hot weather and included significant SNPs from the GWAS. In addition, RHM revealed significant peaks on chromosomes 1 and 13 that are associated with length of productive life and lifetime milk production, respectively, but where GWAS had not identified individual significant SNPs.

Our genomic analyses revealed regions accounting for a detectable variation in the studied traits (h2 = 0.05–0.07 for resilience under hot weather, h2 = 0.05 for lifetime milk production and h2 = 0.06 for length of productive life). Nevertheless, the majority of the trait variance was always accounted for by the polygenic effect. These results are consistent with previous studies, which suggest that fitness-related traits, such as resilience, are usually complex polygenic traits [24]. In the present study, the polygenic effect explained 15% of the variance of resilience to hot weather and 10 to 30% of the variance of the other traits.

Following the genomic analyses of animal resilience to climate change, we performed a functional enrichment analysis, which identified a set of 10 genes encoding olfactory receptors that are involved in multiple biological processes including detection of chemical stimuli in sensory perception and olfactory receptor activity and transduction (see Additional file 9: Fig. S6). Olfactory receptors are expressed in olfactory sensory neurons and detect volatile odorants in smells. However, some receptors are also expressed in other tissues [75, 76]. Response of small ruminants (sheep and goats) to elevated ambient temperatures is largely directed to respiratory evaporation through the respiratory tract in order to dissipate the increased heat load [71, 77, 78]. Heat stress has been related with inducible hypoxia [79]. The OAR5_38936884.1 SNP that we detected on chromosome 5 is located 157 kb upstream the annotated HIG1 gene, which is a member of the hypoxia-inducible domain family and encodes a protein participating in the mitochondrial respiratory chain that catalyses the reduction of oxygen to water. A previous study [80] reported transcriptomic results on the activation of olfactory receptors during hypoxia, which leads to the stimulation of hyperventilation (increased respiration rate) and consequent respiratory evaporation and heat loss. Previously reported findings [81] proposed that Hsc70t, a variant of the heat shock protein 70 (Hsp70), which plays a key role in mitigating heat stress, enhances the expression of mammalian olfactory receptors. Therefore, we suggest that the above genes may constitute the hypoxia sensor in the breathing circuit that increases the respiration rate of the animals and results in thermal homeostasis through respiratory evaporation. Several other candidate genes were also grouped in annotation enrichment clusters (see Additional file 10: Table S4) that are associated with biological processes linked to animal endurance to heat stress [56]. Further research on the transcription of these genes may shed new light into the molecular mechanism that controls sheep resilience to heat stress.


Our results offer a genome-wide understanding of sheep performance resilience to changing weather as a promising adaptive trait for future inclusion in selective breeding programmes. We have identified several genes that contribute to animal milk performance and overall fitness, and genomic regions that are linked to local adaptation of domestic sheep. We propose the integration of different genomic approaches to unravel the genomic regions that are important for breeding purposes as well as for the conservation of genetic variation within livestock breeds in terms of climate change mitigation and improvement of sustainability.

Availability of data and materials

Data analysed during the current study are available from the corresponding author on reasonable request. Phenotypic and genotypic datasets generated are deposited in Mendeley public repository available at


  1. Eusebi PG, Martinez A, Cortes O. Genomic tools for effective conservation of livestock breed diversity. Diversity. 2020;12:8.

    Article  Google Scholar 

  2. Sgrò CM, Lowe AJ, Hoffmann AA. Building evolutionary resilience for conserving biodiversity under climate change. Evol Appl. 2011;4:326–37.

    Article  PubMed  Google Scholar 

  3. FAO. Coping with climate change - the roles of genetic resources for food and agriculture. Rome: FAO; 2015.

  4. Bruford MW, Ginja C, Hoffmann I, Joost S, Orozco-terWengel P, Alberto FJ, et al. Prospects and challenges for the conservation of farm animal genomic resources, 2015–2025. Front Genet. 2015;6:314.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Gunderson KL, Steemers FJ, Lee G, Mendoza LG, Chee MS. A genome-wide scalable SNP genotyping assay using microarray technology. Nat Genet. 2005;37:549–54.

    Article  CAS  PubMed  Google Scholar 

  6. Matukumalli LK, Lawley CT, Schnabel RD, Taylor JF, Allan MF, Heaton MP, et al. Development and characterization of a high density SNP genotyping assay for cattle. PLoS One. 2009;4:e5350.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Bosse M, Megens HJ, Madsen O, Paudel Y, Frantz LA, Schook LB, et al. Regions of homozygosity in the porcine genome: consequence of demography and the recombination landscape. PLoS Genet. 2012;8:e1003100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Scraggs E, Zanella R, Wojtowicz A, Taylor JF, Gaskins CT, Reeves JJ, et al. Estimation of inbreeding and effective population size of full-blood Wagyu cattle registered with the American Wagyu Cattle Association. J Anim Breed Genet. 2014;131:3–10.

    Article  CAS  PubMed  Google Scholar 

  9. 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:71.

    Article  PubMed  PubMed Central  Google Scholar 

  10. 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 (Basel). 2020;10:2285.

    Article  Google Scholar 

  11. 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:220–34.

    Article  CAS  PubMed  Google Scholar 

  12. 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:84.

    Article  PubMed  PubMed Central  Google Scholar 

  13. 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:237–49.

    Article  PubMed  PubMed Central  Google Scholar 

  14. 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:e0176780.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Gorssen W, Meyermans R, Buys N, Janssens S. SNP genotypes reveal breed substructure, selection signatures and highly inbred regions in Pietrain pigs. Anim Genet. 2020;51:32–42.

    Article  CAS  PubMed  Google Scholar 

  16. Gorssen W, Meyermans R, Janssens S, Buys N. A publicly available repository of ROH islands reveals signatures of selection in different livestock and pet species. Genet Sel Evol. 2021;53:2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Williams JL, Hall SJ, Del Corvo M, Ballingall KT, Colli L, Ajmone Marsan P, et al. Inbreeding and purging at the genomic Level: the Chillingham cattle reveal extensive, non-random SNP heterozygosity. Anim Genet. 2016;47:19–27.

    Article  CAS  PubMed  Google Scholar 

  18. Ferenčaković M, Banadinović M, Mercvajler M, Khayat-Zadeh N, Mészáros G, Cubric-Curik V, et al. Mapping of heterozygosity rich regions in Austrian Pinzgauer cattle. Acta Agric Slov. 2016;5:S41–4.

    Google Scholar 

  19. Bizarria W, Pimenta G, Fonseca MG, Pereira GL, Loyola LA, Machado OR, et al. Fine-scale estimation of inbreeding rates, runs of homozygosity and genome-wide heterozygosity levels in the Mangalarga Marchador horse breed. J Anim Breed Genet. 2021;138:161–73.

    Article  Google Scholar 

  20. Joost S, Bruford MW. Editorial: advances in farm animal genomic resources. Front Genet. 2015;6:333.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Berghof TVL, Poppe M, Mulder HA. Opportunities to improve resilience in animal breeding programs. Front Genet. 2019;9:692.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Sánchez-Molano E, Kapsona VV, Oikonomou S, McLaren A, Lambe N, Conington J, et al. Breeding strategies for animal resilience to weather variation in meat sheep. BMC Genet. 2020;21:116.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Macciotta NPP, Biffani S, Bernabucci U, Lacetera N, Vitali A, Ajmone-Marsan P, et al. Derivation and genome-wide association study of a principal component-based measure of heat tolerance in dairy cattle. J Dairy Sci. 2017;100:4683–97.

    Article  CAS  PubMed  Google Scholar 

  24. Sánchez-Molano E, Kapsona VV, Ilska JJ, Desire S, Conington J, Mucha S, et al. Genetic analysis of novel phenotypes for farm animal resilience to weather variability. BMC Genet. 2019;20:84.

    Article  PubMed  PubMed Central  Google Scholar 

  25. International Sheep Genomics Consortium. 2015. Accessed 3 Apr 2015.

  26. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. International Sheep Genomics Consortium, Archibald AL, Cockett NE, Dalrymple BP, Faraut T, Kijas JW, et al. The sheep genome reference sequence: a work in progress. Anim Genet. 2010;41:449–53.

  28. Biscarini F, Cozzi P, Gaspa G, Marras G. detectRUNS: an R package to detect runs of homozygosity and runs of eterozygosity in diploid genomes. R package version 0.9. 5. 2018. Accessed 2 Dec 2020.

  29. Marras G, Gaspa G, Sorbolini S, Dimauro C, Ajmone-Marsan P, Valentini A, et al. Analysis of runs of homozygosity and their relationship with inbreeding in five cattle breeds farmed in Italy. Anim Genet. 2015;46:110–21.

    Article  CAS  PubMed  Google Scholar 

  30. Meyermans R, Gorssen W, Buys N, Janssens S. How to study runs of homozygosity using PLINK? A guide for analyzing medium density SNP data in livestock and pet species. BMC Genomics. 2020;21:94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. 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:359–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies. Nat Genet. 2012;44:821–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Tsartsianidou V, Kapsona VV, Sánchez-Molano E, Basdagianni Z, Carabaño MJ, Chatziplis D, et al. Understanding the seasonality of performance resilience to climate volatility in Mediterranean dairy sheep. Sci Rep. 2021;11:1889.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. 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.

  35. RStudio Team. RStudio: Integrated Development for R. Boston: RStudio, Inc.; 2020. Accessed 20 Aug 2020.

  36. Basdagianni Z, Banos G, Abas Z, Arsenos G, Valergakis GE, Zygoyiannis D. Estimation of daily and total lactation milk yield of Chios ewes from single morning or evening records. Livest Prod Sci. 2005;92:59–68.

    Article  Google Scholar 

  37. Cebamanos L, Gray A, Stewart I, Tenesa A. Regional heritability advanced complex trait analysis for GPU and traditional parallel architectures. Bioinformatics. 2014;30:1177–9.

    Article  CAS  PubMed  Google Scholar 

  38. Gilmour AR, Gogel BJ, Cullis BR, Welham SJ, Thompson R. ASRemlUser guide Release 4.1 Structural specifcation. Hemel Hempstead: VSN International Ltd.; 2015.

  39. Amin N, van Duijn CM, Aulchenko YS. A genomic background based method for association analysis in related individuals. PLoS One. 2007;2:e1274.

  40. Barrett JC, Fry B, Maller J, Daly MJ. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics. 2005;21:263–5.

    Article  CAS  PubMed  Google Scholar 

  41. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GR, Thormann A, et al. The Ensembl Variant Effect Predictor. Genome Biol. 2016;17:122.

    Article  PubMed  PubMed Central  Google Scholar 

  42. The UniProt Consortium. UniProt: the universal protein knowledgebase in 2021. Nucleic Acids Res. 2021;49:D480–9.

    Article  Google Scholar 

  43. Raudvere U, Kolberg L, Kuzmin I, Arak T, Adler P, Peterson H, et al. g:Profiler:a web server for functional enrichement analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 2019;47:W191–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. da Huang W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1–13.

    Article  Google Scholar 

  45. Browning SR, Browning BL. Identity by descent between distant relatives: detection and applications. Annu Rev Genet. 2012;46:617–33.

    Article  CAS  PubMed  Google Scholar 

  46. Zhang Z, Liu Q, Di R, Hu W, Wang X, He X, et al. Single nucleotide polymorphisms in BMP2 and BMP7 and the association with litter size in Small Tail Han sheep. Anim Reprod Sci. 2019;204:183–92.

    Article  CAS  PubMed  Google Scholar 

  47. Dixit SP, Singh S, Ganguly I, Bhatia AK, Sharma A, Kumar NA, et al. Genome-wide runs of homozygosity revealed selection signatures in Bos indicus. Front Genet. 2020;11:92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Suarez-Vega A, Gutierrez-Gil B, Arranz JJ. Transcriptome expression analysis of candidate milk genes affecting cheese-related traits in 2 sheep breeds. J Dairy Sci. 2016;99:6381–90.

    Article  CAS  PubMed  Google Scholar 

  49. Mwacharo JM, Kim ES, Elbeltagy AR, Aboul-Naga AM, Rischkowsky BA, Rothschild MF. Genomic footprints of dryland stress adaptation in Egyptian fat-tail sheep and their divergence from East African and western Asia cohorts. Sci Rep. 2017;7:17647.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Clarkson AH, Paine S, Martin-Tereso J, Kendall NR. Copper physiology in ruminants: trafficking of systemic copper, adaptations to variation in nutritional supply and thiomolybdate challenge. Nutr Res Rev. 2020;33:43–9.

    Article  CAS  PubMed  Google Scholar 

  51. Fraldi A, Biffi A, Lombardi A, Visigalli I, Pepe S, Settembre C, et al. SUMF1 enhances sulfatase activities in vivo in five sulfatase deficiencies. Biochem J. 2007;403:305–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Marina H, Reverter A, Gutiérrez-Gil B, Alexandre PA, Porto-Neto LR, Suarez-Vega A, et al. Gene networks driving genetic variation in milk and cheese-making traits of Spanish Assaf sheep. Genes (Basel). 2020;11:715.

    Article  CAS  Google Scholar 

  53. Zhang F, Qu K, Chen N, Hanif Q, Jia Y, Huang Y, et al. Genome-wide SNPs and InDels characteristics of three Chinese cattle breeds. Animals (Basel). 2019;9:596.

    Article  CAS  Google Scholar 

  54. Li X, Yang J, Shen M, Xie XL, Liu GJ, Xu YX, et al. Whole-genome resequencing of wild and domestic sheep identifies genes associated with morphological and agronomic traits. Nat Commun. 2020;11:2815.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Bernabucci U, Lacetera N, Baumgard LH, Rhoads RP, Ronchi B, Nardone A. Metabolic and hormonal acclimation to heat stress in domesticated ruminants. Animal. 2010;4:1167–83.

    Article  CAS  PubMed  Google Scholar 

  56. Carabaño MJ, Ramón M, Menéndez-Buxadera A, Molina A, Díaz C. Selecting for heat tolerance Anim Front. 2019;9:62–8.

    Article  PubMed  Google Scholar 

  57. Kijas JW, Lenstra JA, Hayes B, Boitard S, Porto Neto LR, San Cristobal M, et al. Genome-wide analysis of the world’s sheep breeds reveals high levels of historic mixture and strong recent selection. PLoS Biol. 2012;10:e1001258.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Mastrangelo S, Di Gerlando R, Tolone M, Tortorici L, Sardina MT, Portolano B, et al. Genome wide linkage disequilibrium and genetic structure in Sicilian dairy sheep breeds. BMC Genet. 2014;15:108.

    Article  PubMed  PubMed Central  Google Scholar 

  59. Antonios S, Rodríguez-Ramilo ST, Aguilar I, Astruc JM, Legarra A, Vitezica ZG. Genomic and pedigree estimation of inbreeding depression for semen traits in the Basco-Bearnaise dairy sheep breed. J Dairy Sci. 2021;104:3221–30.

    Article  CAS  PubMed  Google Scholar 

  60. Meyermans R, Gorssen W, Wijnrocx K, Lenstra JA, Vellema P, Buys N, et al. Unraveling the genetic diversity of Belgian Milk Sheep using medium-density SNP genotypes. Anim Genet. 2020;51:258–65.

    Article  CAS  PubMed  Google Scholar 

  61. Pekkala N, Knott KE, Kotiaho JS, Nissinen K, Puurtinen M. The effect of inbreeding rate on fitness, inbreeding depression and heterosis over a range of inbreeding coefficients. Evol Appl. 2014;7:1107–19.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Forutan M, Ansari Mahyari S, Baes C, Melzer N, Schenkel FS, Sargolzaei M. Inbreeding and runs of homozygosity before and after genomic selection in North American Holstein cattle. BMC Genomics. 2018;19:98.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Gutiérrez-Gil B, Arranz JJ, Pong-Wong R, García-Gámez E, Kijas J, Wiener P. Application of selection mapping to identify genomic regions associated with dairy production in sheep. PLoS One. 2014;9:e94623.

    Article  PubMed  PubMed Central  Google Scholar 

  64. Curik I, Ferenčaković M, Sölkner J. Inbreeding and runs of homozygosity: A possible solution to an old problem. Livest Sci. 2014;166:26–34.

    Article  Google Scholar 

  65. Kim ES, Elbeltagy AR, Aboul-Naga AM, Rischkowsky B, Sayre B, Mwacharo JM, et al. Multiple genomic signatures of selection in goats and sheep indigenous to a hot arid environment. Heredity (Edinb). 2016;116:255–64.

    Article  CAS  Google Scholar 

  66. Yurchenko AA, Deniskova TE, Yudin NS, Dotsev AV, Khamiruev TN, Selionova MI, et al. High-density genotyping reveals signatures of selection related to acclimation and economically important traits in 15 local sheep breeds from Russia. BMC Genomics. 2019;20:S294.

    Article  Google Scholar 

  67. Guan D, Luo N, Tan X, Zhao Z, Huang Y, Na R, et al. Scanning of selection signature provides a glimpse into important economic traits in goats (Capra hircus). Sci Rep. 2016;6:36372.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Thorgeirsdottir S, Sigurdarson S, Thorisson HM, Georgsson G, Palsdottir A. PrP gene polymorphism and natural scrapie in Icelandic sheep. J Gen Virol. 1999;80:2527–34.

    Article  CAS  PubMed  Google Scholar 

  69. Goldmann W, Houston F, Stewart P, Perucchini M, Foster J, Hunter N. Ovine prion protein variant A136R154L168Q171 increases resistance to experimental challenge with bovine spongiform encephalopathy agent. J Gen Virol. 2006;87:3741–5.

    Article  CAS  PubMed  Google Scholar 

  70. Finocchiaro R, van Kaam JB, Portolano B, Misztal I. Effect of heat stress on production of Mediterranean dairy sheep. J Dairy Sci. 2005;88:1855–64.

    Article  CAS  PubMed  Google Scholar 

  71. Sevi A, Caroprese M. Impact of heat stress on milk production, immunity and udder health in sheep: A critical review. Small Ruminant Res. 2012;107:1–7.

    Article  Google Scholar 

  72. Giorgi F, Bi X. Updated regional precipitation and temperature changes for the 21st century from ensembles of recent AOGCM simulations. Geophys Res Lett. 2005;32:L21715.

    Article  Google Scholar 

  73. Carabaño MJ, Bachagha K, Ramón M, Díaz C. Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools. J Dairy Sci. 2014;97:7889–904.

    Article  PubMed  Google Scholar 

  74. Nagamine Y, Pong-Wong R, Navarro P, Vitart V, Hayward C, Rudan I, et al. Localising loci underlying complex trait variation using Regional Genomic Relationship Mapping. PLoS One. 2012;7:e46501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Fleischer J, Breer H, Strotmann J. Mammalian olfactory receptors. Front Cell Neurosci. 2009;3:9.

    Article  PubMed  PubMed Central  Google Scholar 

  76. Kang N, Koo J. Olfactory receptors in non-chemosensory tissues. BMB Rep. 2012;45:612–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Hamzaoui S, Salama AA, Albanell E, Such X, Caja G. Physiological responses and lactational performances of late-lactation dairy goats under heat stress conditions. J Dairy Sci. 2013;96:6355–65.

    Article  CAS  PubMed  Google Scholar 

  78. Marai IFM, El-Darawany AA, Fadiel A, Abdel-Hafez MAM. Physiological traits as affected by heat stress in sheep—A review. Small Ruminant Res. 2007;71:1–12.

    Article  Google Scholar 

  79. Paul C, Teng S, Saunders PT. A single, mild, transient scrotal heat stress causes hypoxia and oxidative stress in mouse testes, which induces germ cell death. Biol Reprod. 2009;80:913–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Chang AJ, Ortega FE, Riegler J, Madison DV, Krasnow MA. Oxygen regulation of breathing through an olfactory receptor activated by lactate. Nature. 2015;527:240–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Neuhaus EM, Mashukova A, Zhang W, Barbour J, Hatt H. A specific heat shock protein enhances the expression of mammalian olfactory receptor proteins. Chem Senses. 2006;31:445–52.

    Article  CAS  PubMed  Google Scholar 

Download references


We thank Chios Sheep Breeder’s Cooperative Macedonia for animal records provided. We also thank the Hellenic National Meteorological Service for the availability of weather data. We acknowledge the contribution of Dr. Sotiria Vouraki and Angeliki Argyriadou, School of Veterinary Medicine, Aristotle University of Thessaloniki, to blood sample collection.


The work was mainly funded by iSAGE Horizon 2020 project: “Innovation for Sustainable Sheep and Goat Production in Europe” (679302; Additional support was received by the UK Biotechnology and Biological Sciences Research Council ISP3 (BBS/E/D/30002275) and the Rural & Environment Science & Analytical Services (RESAS) Division of the Scottish Government.

Author information

Authors and Affiliations



ZB organized the animal phenotypic data using the database of Chios Sheep Breeder’s Cooperative Macedonia. GA and DC provided blood and DNA samples, respectively. VT, GB and ZB performed the animal data quality control. All analyses were performed by VT supported by VVK and ESM, under the supervision of GB and DC. VT visualised data and results and drafted the first manuscript, and GB reviewed the first draft. All authors reviewed and edited the ensuing and final draft. GB, AT and GA supervised the project and acquired funding. GB was the principal investigator. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Valentina Tsartsianidou.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Correlation heatmap between different measures of inbreeding coefficient. The respective coefficients are: SNP-based genomic inbreeding (F), whole-genome ROH-based inbreeding (FROH), and FROH by the respective size classes (ROH segment length) 2 to 4 Mb (FROH2_4), 4 to 8 Mb (FROH4_8), 8 to 16 Mb (FROH8_16), >16 Mb (FROH16).

Additional file 2: Table S1.

Islands of runs of homozygosity (ROH) and heterozygosity (ROHet), and associated genes and animal traits of the studied sheep. The respective chromosomal locations are attributed to Oar v4.0 genome assembly.

Additional file 3

: Table S2. Functional characterization of annotated genes in detected islands of runs of homozygosity (ROH) and heterozygosity (ROHet).

Additional file 4: Figure S2.

Examples of individual animal resilience phenotypes manifested in reaction norm curves. Each curve represents the individual animal change in daily milk yield (milk, kg) in response to average air temperature variation on the day of milk record (Temperature, oC).

Additional file 5: Figure S3.

Quantile-Quantile plots for genome-wide association. Figures show the Q-Q plot before correction by the inflation factor for animal resilience to temperature change in hot weather (a), lifetime milk yield (b), length of productive life (c) and cold weather (d).

Additional file 6

: Figure S4. Linkage disequilibrium (r2) heatmap of chromosome 5 region associated with Chios sheep resilience to hot weather conditions. This region extends between the furthest SNPs above the genome-wide and suggestive significance thresholds spanning a region of 2.41 Mb for animal resilience to temperature change in hot weather. The pairwise LD between these two SNPs is illustrated at the bottom right part of the figure.

Additional file 7: Table S3.

List of the genes identified in the genomic analysis of resilience to temperature change under hot weather conditions. The respective genes were detected 1 Mb upstream and downstream of the genome-wide and suggestive significant SNPs.

Additional file 8:

Figure S5. Protein interaction network. The present network represents the candidate identified genes for sheep resilience to temperature fluctuations under hot weather.

Additional file 9:

Figure S6. g:GOSt multi-query Manhattan plot of functional enrichment analysis for sheep resilience to temperature fluctuations under hot weather. The X-axis represents the functional terms colour-separated by data source. The inset table summarizes the biological processes of the genes.

Additional file 10:

Table S4. Significantly enriched annotation clusters and functional terms.

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 The Creative Commons Public Domain Dedication waiver ( 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

Tsartsianidou, V., Sánchez-Molano, E., Kapsona, V.V. et al. A comprehensive genome-wide scan detects genomic regions related to local adaptation and climate resilience in Mediterranean domestic sheep. Genet Sel Evol 53, 90 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: