Skip to main content

Using high-density SNP data to unravel the origin of the Franches-Montagnes horse breed

Abstract

Background

The Franches-Montagnes (FM) is the last native horse breed of Switzerland, established at the end of the 19th century by cross-breeding local mares with Anglo-Norman stallions. We collected high-density SNP genotype data (Axiom™ 670 K Equine genotyping array) from 522 FM horses, including 44 old-type horses (OF), 514 European Warmblood horses (WB) from Sweden and Switzerland (including a stallion used for cross-breeding in 1990), 136 purebred Arabians (AR), 32 Shagya Arabians (SA), and 64 Thoroughbred (TB) horses, as introgressed WB stallions showed TB origin in their pedigrees. The aim of the study was to ascertain fine-scale population structures of the FM breed, including estimation of individual admixture levels and genomic inbreeding (FROH) by means of Runs of Homozygosity.

Results

To assess fine-scale population structures within the FM breed, we applied a three-step approach, which combined admixture, genetic contribution, and FROH of individuals into a high-resolution network visualization. Based on this approach, we were able to demonstrate that population substructures, as detected by model-based clustering, can be either associated with a different genetic origin or with the progeny of most influential sires. Within the FM breed, admixed horses explained most of the genetic variance of the current breeding population, while OF horses only accounted for a small proportion of the variance. Furthermore, we illustrated that FM horses showed high TB admixture levels and we identified inconsistencies in the origin of FM horses descending from the Arabian stallion Doktryner. With the exception of WB, FM horses were less inbred compared to the other breeds. However, the relatively few but long ROH segments suggested diversity loss in both FM subpopulations. Genes located in FM- and OF-specific ROH islands had known functions involved in conformation and behaviour, two traits that are highly valued by breeders.

Conclusions

The FM remains the last native Swiss breed, clearly distinguishable from other historically introgressed breeds, but it suffered bottlenecks due to intensive selection of stallions, restrictive mating choices based on arbitrary definitions of pure breeding, and selection of rare coat colours. To preserve the genetic diversity of FM horses, future conservation managements strategies should involve a well-balanced selection of stallions (e.g., by integrating OF stallions in the FM breeding population) and avoid selection for rare coat colours.

Background

The Franches-Montagnes (FM) is the last native horse breed of Switzerland, established at the end of the 19th century from crossbreeding local Jura and Comtois mares with two Anglo-Norman founding sires, Vaillant (1891) and Imprévu (1886). The FM studbook was officially released in 1921 [1, 2]. Periodic introgression of horses from other breeds were allowed to meet market demands, initially for heavier draught horses required for agricultural work, which led to cross-breeding with Percheron, Breton, and Belgian draught horses [1, 2]. In response to the increased mechanisation of agriculture after World War II, additional introgressions using Shagya Arabian, Purebred Arabian, Nonius, Anglo-Norman, Selle Français, and Swiss Warmblood stallions were carried out to select a lighter horse type for leisure riding and carriage driving. Some of these introgressions also led to the formation of new sire lines. Traditionally, the sire lines of the FM breed are expressed as letters, where each stallion’s name starts with the same letter as his  sire. Currently, there are 11 extant sire lines: C, D, E, H, R, and V are considered original, although the most influential V sire (Vulcain, 1983) had a Warmblood damsire. Additional lines, such as the DON, L, N, P, and Q lines, descend from sires used for cross-breeding. The DON line originated from the purebred Arabian stallion Doktryner (1950), which covered FM mares from 1958, with only one of his sons (Don Pablo, 1964) producing licenced stallions himself. In the 1970s, the Swedish Warmblood stallions Aladin (1964) and Nello (1971), sons of the Trakhener stallion Nepal (1956) were selected, resulting in the sire lines L and P. The last introgressions involved two Swiss Warmblood stallions from the same dam (Salomé, 1977), Qui-Sait (1985) and Noé (1984), establishing the sire lines Q and N, respectively.

To preserve the breed characteristics obtained during its admixed history, the studbook was closed in 1997. To date, admixture proportions are estimated based on pedigree information as a total percentage, while breed origin is ignored. Furthermore, by convention, FM horses born before 1950 are considered purebred, meaning that earlier introgressions are not accounted for in the current computation of admixture proportions. Since 2012, purebred or old-type FM (hereafter, OF and considered as a separate breed) are managed as a subsection of the FM studbook. The OF selection criteria (field test, stallion licensing) differ from those of the main FM population. The main breeding goal of the OF is to maintain a heavier draught horse, while the focus of the FM is to breed a leisure horse for riding and carriage driving. Both FM and OF breeders share a common goal of preserving the perceived good temperament of the breed. Currently, OF is declining in numbers, resulting in increased inbreeding in this subpopulation. Conversely, additional cross-breeding activities threaten the integrity of the FM horse population as an original native breed, both from a genetic [3] and a phenotypic perspective [4].

Runs of homozygosity (ROH) in an individual’s genome, which consist of haplotypes identical-by-descent (IBD), are one of the main tools to quantify genomic levels of inbreeding. The length of ROH segments can also be used to ascertain historical changes in population size and structure. Recombination and admixture events are major reasons for breakdown of ROH, such that long ROH segments indicate recent inbreeding, while the genome-wide assessment of ROH segments has emerged as a valuable tool to estimate inbreeding. Hence, outcrossed animals typically carry fewer and shorter ROH segments, while animals descending from long-term inbred populations carry more and shorter ROH segments. Recently inbred animals (e.g. bottleneck) carry fewer and longer ROH segments (see Ceballos et al. [5] for a complete review). The genomic inbreeding coefficient (FROH) for an animal is derived by dividing the sum of all homozygous segments (SROH) by the total length of the genome [6]. Numerous studies have demonstrated that same overlapping ROH segments between multiple individuals, so-called ROH islands, also highlight genetic regions that have been under selection, in cattle [7, 8], sheep [9,10,11], and honeybees [12]. In horses, based on different genotyping platforms and ROH settings, numerous breed-specific ROH islands have already been identified that include genes associated with behaviour [13, 14], coat colour [14,15,16,17,18], coat quality [19], fertility [13, 18,19,20,21,22], size [15, 22], skeletal development [18, 19, 23], muscle growth [24], geographical adaptations [14, 25], and gait quality [14].

In this study, we used high-density SNP genotypes from the Affymetrix Axiom Equine genotyping array to ascertain the fine-scale population structure of the FM breed, taking into account the last four cross-breeding activities. Furthermore, we determined the level of FROH in all breeds by means of ROH. Within FM we additionally compared FROH against pedigree-based inbreeding (FPED) and admixture proportions. Finally, we identified private and overlapping ROH islands to determine genomic regions that are under selection in the respective breed.

Methods

Dataset

The genomic data consisted of 1268 horses originating from five different European breeds. We collected 44 samples of OF and 478 of FM. Between 2014 and 2023, we genotyped nearly all three-year old males selected for the stallion licensing test, since genetic testing for congenital hepatic fibrosis became mandatory in 2012 [26]. The DNA from older stallions was available through frozen sperm, including  two F1 sires that descended from the Warmblood stallions Noé and Qui-Sait. Additional geldings and mares were genotyped in the same period of time, such that all 11 extant sire lines were represented by eight (Q line) up to 148 descendants (H line).

Furthermore, we included 32 Shagya Arabians (SA) because of the Shagya introgression in the 1940s, 136 Purebred Arabians (AR) descending from the introgressed stallion Doktryner at the 5th to 8th generation, and 514 European Warmblood horses (WB, 380 belonging to the Swedish Warmblood Studbook, and 134 to the Swiss Warmblood Studbook). The WB samples included one of the two Swiss Warmblood stallions (Noé, 1984) that was used for cross-breeding in the 1990s. Finally, we included a convenience sample of 64 Thoroughbred horses (TB) from Australasia, as the WB stallions used for cross-breeding also showed TB origin in their pedigrees.

The DNA from FM and Swiss WB was sampled under permits VD3527b, VD2976.1, and VD2227.2, each approved by the cantonal veterinary office of Vaud, Switzerland, while the genotype information from the other breeds was derived from previous studies [13, 14].

SNP genotyping

All 1268 horses were genotyped on the commercial Axiom™ Equine Genotyping array that contained 670,796 evenly distributed single nucleotide polymorphisms (SNPs) [27]. The chromosomal position of the SNPs was determined based on the EquCab3.0 reference genome [28] and the most recent annotation file (https://www.thermofisher.com/order/catalog/product/550583). Single nucleotide polymorphisms positioned on the sex chromosomes or with unknown chromosomal position were excluded from further analyses, resulting in 602,131 SNPs. Missing genotype information was imputed separately for each breed using BEAGLE 5.2 [29]. For population clustering and admixture analysis, we excluded SNPs with a minor allele frequency (MAF) less than 0.05, resulting in 421,988 SNPs per horse.

Population clustering

To illustrate the population structure of the horses, we performed a principal component analysis (PCA) and computed pairwise FST distances between the breeds. Principal Component Analysis was applied to a genetic relationship matrix (G) with pairwise identities by state (IBS) between horses, as provided by PLINK v1.9 [30], and the first two principal components (PCs) were visualized using the R-platform (www.r-project.org) [31]. We also calculated pairwise FST estimates between breeds using the Wright’s FST function, as implemented in PLINK v 1.9 [30], and used the program SPLITSTREE (www.splistree.org) [32] to create phylogenetic networks between the breeds.

Admixture

Admixture proportions for each horse were determined using the program Admixture 1.23 [33]. We ran Admixture for 100 iterations, increasing K from 2 to 20. Convergence between independent runs at the same K was monitored by comparing the resulting log-likelihood scores (LLs) following 100 iterations and was inferred from stabilized LLs with less than 1 LL unit of difference between runs. Cross validation (CV) error estimation for each K was performed to determine the optimal number of clusters. We also derived the admixture proportions of FM horses using the supervised learning mode, by specifying the ancestries of the reference individuals according to the five sampled breeds (TB, AR, SA, WB and OF). Supervised (K = 5) and unsupervised admixture results by increasing K from 2 to 10, were visualized with the program Distruct 1.1 [34] and integrated in the high-resolution FM population network, as described below.

Identification of key contributors

Within the FM breed, in addition to the admixture proportions, we also identified key contributors to ascertain the impact of admixed horses on the current population. The method to identify key contributors which explain most of the genetic variance in complex population structures was described in detail by Neuditschko et al. [35] and available online at https://github.com/esteinig/netview. It requires as input a symmetric relationship matrix A of dimension n × n and the number of k significant PCs, where n is the number of individuals. Here, we used the aforementioned IBS relationship matrix G and determined the number of k significant PCs with the empirical method described as Horn’s parallel analysis, which is implemented in the statistical software package paran [36]. This method employs Monte Carlo estimates to retain the most significant PCs under a defined level of significance and number of iterations. Here, we chose a significance level of P = 0.01 and 10,000 iterations, which have been suggested in the modified version of Horn’s parallel analysis [37]. After determining the number of k significant PCs, we calculated the genetic contribution score (gcj) for each horse. The computation of the gcj is based on a Singular Value Decomposition (SVD) of G and accounts for the correlation between the jth individual relationship vector and the ith standardized eigenvector, limiting the number of eigenvectors to the first k significant PCs. Thereafter, we ranked the horses according to gcj and also explored the association between gcj and individual admixture proportions.

High-resolution population networks

A high-resolution population network visualization of FM horses was performed by applying NetView, as described in detail by Neuditschko et al. [38] and Steinig et al. [39]. Briefly, we computed genetic distances by subtracting pairwise IBS relationships from 1 and applied the algorithm in its default setting (number of k nearest neighbours k-NN = 10). To illustrate the genetic relatedness between neighbouring horses, we associated the thickness of edges (connecting lines) with the magnitude of the genetic distance, with thicker edges corresponding to lower genetic distances. To identify highly influential and inbred horses within the respective population networks, we scaled the node size of each horse based on the individual gcj and FROH (see computation below), with the node colour representing the individual admixture proportions at the respective number of clusters (unsupervised K = 7 and supervised K = 5).

Runs of homozygosity

Runs of homozygosity segments were determined with an overlapping window approach implemented in PLINK v1.9 [30], using all 602,131 genome-wide mapped SNPs. The following settings were applied: a minimum SNP density of one SNP per 50 kb, a maximum gap length of 100 kb, and a minimum length of homozygous segments of 500 kb (including more than 80 homozygous SNP genotypes), while one heterozygous SNP genotype was permitted per segment. The total number of ROH (NROH), the total length of ROH segments (SROH), and the average length of ROH (LROH) were summarised for each breed. The genomic inbreeding coefficients (FROH) were calculated by dividing SROH by the length of the autosomal genome (LAUTO), which was set to 2280.92 Mb. Differences between breeds were investigated using ANOVA and post hoc Tukey’s HSD (honestly significant difference) tests at a significance level of α < 0.05, as implemented in the R package multcompView [40]. The length of segments was also categorised for each breed to explore breed-specific ROH islands. We also compared FROH and cumulative genomic admixture proportions (AGEN) of each FM horse with their respective pedigree-based inbreeding (FPED) and admixture proportions (APED), using a linear regression model, as implemented in the statistical computing software R [31]. The cumulative genomic admixture proportions were provided by the FM association, while FPED of the horses was calculated with the R package pedigree [41].

Runs of homozygosity islands and gene functions

Putative ROH islands were determined based on overlapping homozygous regions that were shared by more than 50% of horses within each breed, using the R package detectRUNS [42]. Runs of homozygosity islands that occurred in only one breed were defined as “breed-specific” or “private”. Additionally, we also determined overlapping ROH islands between breeds, using a threshold for the minimal overlap of at least two SNPs. We used the NCBI genome data viewer (https://www.ncbi.nlm.nih.gov/genome/gdv/) and the reference genome assembly EquCab3.0 [28] to identify genes located in ROH islands. Gene ontology (GO) analysis on the identified genes was performed for each breed using the term enrichment service provided on https://amigo.geneontology.org/amigo, with the settings reference list “Equus caballus” and Fisher’s Exact test with Bonferroni correction for multiple testing. The gene list was limited to uniquely mapped gene IDs. Known functions of the identified genes were assigned  by conducting a literature review and allele frequencies of functional SNPs were calculated for each breed, i.e. SNPs with known effects on specific phenotypes.

Results

Population clustering

Visualization of the first versus the second PC showed that the six breeds comprise three distinct population clusters (Fig. 1a). The first PC, accounting for 50% of the genetic variance, differentiated OF and FM from the other breeds, while the second PC, accounting for 10% of the genetic variance, separated the two Arabian breeds (AR and SA) from WB and TB. Based on the first two PCs, it was not possible to clearly distinguish the two Arabian breeds from each other, i.e. WB from TB and OF from FM horses. The two F1 crosses showed a noticeable distance to the remaining FM horses, while the WB sire that was used for cross-breeding was not positioned next to the F1 crosses. The topology of the phylogenetic network based on FST distances coincided with the visualization of the first two PCs, thereby refining the genetic distances between breeds (Fig. 1b). The lowest FST distance between breeds was observed for WB and TB (0.05), followed by SA and AR (0.06). Old-type FM horses showed the highest FST distances to all breeds, in the following decreasing order TB (0.19), SA (0.16), AR (0.14), WB (0.11). Compared to OF, pairwise FST distances between FM and the other breeds were slightly lower. The largest differences between OF and FM pairwise FST distances were observed for TB (0.04), followed by SA (0.03), while the distance between the two subpopulations, FM and OF, was 0.02.

Fig. 1
figure 1

Population clustering of the six horse breeds. a Visualization of horses on the first two principal components (PCs) explaining 60% of the genetic variance; and b phylogenetic network of pairwise FST distances between breeds and subpopulations

Admixture

While increasing K from 1 to 20, the CV error estimation did not provide an optimal cluster solution (see Additional file 1: Figure S1). The first level of model-based clustering (K = 2) confirmed the findings of the first PC, by allocating FM and OF horses into a distinct population cluster and simultaneously pinpointed admixed FM horses, while OF horses were less admixed (Fig. 2). At K = 3, the two Arabian breeds (AR and SA) formed a distinct cluster, whereas some AR horses showed a high level of admixture with TB horses. At the next level (K = 4), WB horses were differentiated from TB horses, although all WB horses showed high TB admixture levels. At two additional level of clustering, FM (K = 5) and WB (K = 6), horses were sub-structured, before SA were differentiated from AR horses at K = 7. Hence, we considered K = 7 to be an optimal cluster solution, as all sampled breeds built a distinct population cluster. Further, when increasing K to 10, in particular WB and FM horses were further sub-structured, while at K = 8 and K = 10, WB horses from Switzerland were allocated to a distinct population cluster (light blue), although these horses already showed a different admixture pattern at all previous levels of clustering (e.g., at K = 4 WB horses from Switzerland had higher TB and Arabian admixture levels compared to WB horses from Sweden).

Fig. 2
figure 2

Model-based clustering result. Graphical representation of individual cluster membership coefficients of horses increasing the number of clusters (K) from 2 to 10

Re-ordering the cluster membership proportions of FM horses at K = 7 according to the 11 sire lines (Fig. 3a) also revealed that the FM substructures were associated with a different genetic background, as especially horses that descended from the introgressed stallions Noé (sire line N) and Aladin (sire line L) were assigned to distinct population clusters, while the substructure associated with sire line N (last introgression) was already apparent at K = 5 (Fig. 2). However, it should be noticed, that the admixture profile of FM horses at K = 7 was poorly associated with sire line origin, as horses that originated from a different sire line showed the same admixture pattern and vice versa. Applying the supervised learning mode at K = 5, admixed FM horses were identified across all sire lines, while especially horses from sire lines DON, L, N, and Q showed high TB and WB admixture levels (Fig. 3b). In this context, it was not expected that horses descending from the Arabian stallion Doktryner (line DON) had higher TB than AR admixture levels, which were similar to those observed for horses descending from the Warmblood stallion Noé (line N). All admixed FM horses displayed unanticipated high TB admixture levels, as no TB stallions were used for cross-breeding in the FM breeding history.

Fig. 3
figure 3

Admixture results and high-resolution population network visualizations of FM horses. a Membership coefficients of FM horses, ordered along stallion line, at unsupervised clustering K = 7; b membership coefficients of FM horses, ordered along stallion line, at supervised clustering K = 5; c high resolution population network, where each FM horse is illustrated by a node, with individual node size proportional to \({gc}_{j}\), while the node color represents the individual levels of admixture at unsupervised K = 7. The thickness of edges varies in the proportion of the genetic distance to visualize individual relationships between the horses. Top key contributors are highlighted by their respective stallion line; and d high resolution population network, where each FM horse is illustrated by a node, with individual node size proportional to FROH, while the node colour represents the individual levels of admixture at supervised K = 7. The two F1 sires and the highly inbred grey horse are highlighted by images and their respective stallion line

Identification of key contributors and high-resolution network visualization

Applying Horn’s parallel analysis on \(G\) including only FM horses resulted in 33 significant PCs that accounted for 78% of the genetic variance (result not shown). After determining \({gc}_{j}\) and FROH for each horse, we performed two high-resolution network visualizations (Fig. 3c, d). The network visualization that illustrates \({gc}_{j}\) and individual admixture levels of each horse at unsupervised clustering K = 7, represented a common observed topology, where key contributors and their respective progeny built distinct population clusters at the edge of the network (Fig. 3c). Simultaneously, it was noted that some of these population clusters caused the aforementioned substructures in the FM population. Ranking the horses according to \({gc}_{j}\) showed that the F1 sire Népal was amongst top key contributors (ranked on 11th position), while the other F1 sire Queens (ranked on the 387th position) and especially OF horses, which were positioned at the middle-left of the network, were less influential. In contrast, the population network associated with FROH and individual admixture levels of each horse at supervised clustering K = 5 revealed that OF horses are highly inbred, while in general admixed horses were less inbred (Fig. 3d). However, in this context, a highly inbred grey horse with an Arabian admixture proportion of 3% stood out.

Comparison between pedigree-based and genomic derived inbreeding and admixture proportions of FM horses

The correlation between genomic admixture levels derived from the unsupervised model (K = 2) and from the supervised model (K = 5) was 0.99. Supervised admixture levels (AGEN) were compared with pedigree-based admixture (APED) in FM only, as OF was used as reference. The mean APED was 11.50 ± 6.59, with a minimum of 0.39 and a maximum of 50.00 for the two F1 crosses. The mean AGEN was 6.40 ± 6.26, with a minimum of 0.004 and a maximum of 41.81, nearly half of the mean APED. Despite this large discrepancy, the Pearson’s product-moment correlation between AGEN and APED was 0.82 (R2 = 0.68, Fig. 4a) and highly significant (p < 0.0001). The intercept was − 2.60, meaning that AGEN was generally only 2.60% lower than APED, suggesting individuals with large divergences between the two values. Among the individuals with a higher AGEN than APED, horses with Don Fernando (F2 to Doktryner) ancestry stood out: the largest difference was seen in Don Fernando’s direct son (F3), with 15% greater AGEN than APED.

Fig. 4
figure 4

Comparison between genomic and pedigree-based estimates. a Linear regression of the relationship between the pedigree-based admixture (APED) and the genomic admixture (AGEN) from the supervised model at K = 5; and b linear regression of the relationship between the pedigree-based inbreeding coefficient (FPED) and the genomic inbreeding coefficient FROH. The OF horses are represented in dark green, the FM in lighter green

The mean pedigree-based inbreeding (FPED) was 6.84 ± 1.83 for OF and FM combined, while the mean FROH was 12.15 ± 2.68. The Pearson’s product-moment correlation between FROH and FPED was 0.56 (R2 = 0.31) and highly significant (p < 0.0001). The intercept for the linear regression between FROH and FPED was 6.53, meaning that FROH was generally 6.53% higher than FPED, i.e. nearly double (Fig. 4b). However, by adding the information of sire and dam to the regression equation (FROH ~ FPED + sire + dam), R2 increased to 0.79.

ROH patterns and distribution

The number of ROH segments was significantly different between all breeds, except between OF and WB (Table 1). Franches-Montagnes had the lowest mean number of ROH segments (NROH), TB the highest. Mean total length of ROH segments (SROH) significantly differed between all breeds, except between OF and SA. The shortest mean average length was for TB, the longest for OF. The mean LROH for FM and OF significantly differed from that of the other four breeds. The mean FROH differed between all breeds: TB had the highest mean FROH, followed by AR, SA, OF, FM, and WB. Old-type FM had more ROH segments, but they were not significantly longer on average (LROH) compared to FM. Therefore, OF had a longer SROH and were more inbred than FM.

Table 1 Summary statistics for number (NROH), total length (SROH), and average length (LROH) of ROH segments, and of genomic inbreeding (FROH)

While the total mean and total length of ROH segments was highest in TB, AR had highest proportion of ROH segments with a length > 10 MB, followed by FM and OF (Fig. 5). Warmblood and SA also had a high frequency of short ROH segments (0.5–2 MB), but a slightly higher proportion of ROH segments of intermediate length (4–8 MB) for WB and of long ROH segments (8–10 MB) for SA. Old-type FM had a majority of short to medium length ROH segments, while FM had proportionally more very long segments.

Fig. 5
figure 5

Distribution of the length of runs of homozygosity (ROH) segments. Length classes: 0.5–1, > 1–2, > 2–4, > 4–6, > 6–8, > 8–10 and > 10 Mb. Breeds: TB: Thoroughbred; AR: Purebred Arabian; SA: Shagya Arabian; WB: Warmblood; OF: old-type Franches-Montagnes; FM: Franches-Montagnes

ROH islands

There were no ROH islands on ECA13, ECA20, and ECA24 to ECA31. Thoroughbred had the most ROH islands (n = 78, Table 2), which concurs with results based on NROH, i.e. breeds with more NROH also had more ROH islands. Warmblood had no ROH islands. Old-type FM had the longest ROH island (2.46 Mb). The chromosomal location and genes contained in the ROH islands were summarised by breed (see Additional file 2: Table S1, Additional file 3: Table S2, Additional file 4: Table S3, Additional file 5: Table S4, Additional file 6: Table S5). Not all loci in the ROH islands were annotated. The results from the GO analysis for each breed are summarized in Additional file 7: Table S6, Additional file 8: Table S7, Additional file 9: Table S8, Additional file 10: Table S9, and Additional file 11: Table S10. Most notably, genes involved in the biological process epithelial cell differentiation (GO:0030855) were overrepresented in AR (see Additional file 7: Table S6) and those in the biological process embryonic skeletal system morphogenesis (GO:0048704) in FM (see Additional file 11: Table S10) and OF (see Additional file 10: Table S9).

Table 2 Summary of runs of homozygosity (ROH) islands per breed

There were 24 ROH islands that overlapped between at least two breeds (Table 3). One ROH island on ECA3 was shared between FM and Arabian breeds. Old-type FM shared ROH islands with FM, AR and SA and one with TB; FM shared islands with AR and SA, and the two Arabian breeds shared islands with TB.

Table 3 Overlaps of runs of homozygosity (ROH) islands (beginning and end in bp, length in Mb) between breeds, including annotated genes

Genotype frequencies of functional SNPs

We calculated the allele frequencies of functional SNPs in the following genes: LASP1, UVSSA, MSTN, MC1R, and NTM (Fig. 6). For LASP1, which has previously been associated with body size [15, 43], FM had the highest frequency of the “big” allele, even higher than OF, followed by TB and WB. Warmblood horses had the highest frequency of the UVSSA allele that was previously associated with lower UV exposure and therefore adaptation to Northern climates [14]. Old-type FM horses had the highest frequency of the allele associated with increased UV exposure, e.g. in Mountain ranges. The homozygous state for MC1R, expressing the chestnut coat colour, had the highest frequency AR, OF, FM, and WB. The majority of TB horses were homozygous for the MSTN allele, which has been shown to be beneficial for short distance races (C–C genotype) [44] while all SA and the majority of AR horses were homozygous for the opposite allele. For the intronic SNP for NTM, which was previously associated with the number of race course starts and therefore with precocity in TB [45], OF had the highest proportion of the favourable C allele, while TB had only few homozygous C–C horses. Franches-Montagnes horses with C–C genotypes had on average lower APED compared to C–T and T–T horses (see Additional file 12: Figure S2).

Fig. 6
figure 6

Frequencies of SNP genotypes with known functions in the six breeds  Thoroughbred (TB), Purebred Arabian (AR), Shagya Arabian (SA), European Warmblood (WB), Old-type Franches-Montagnes (OF) and Franches-Montagnes (FM). a The A allele of the LASP1 SNP was associated with large sized horses; b the A allele from the UVSSA SNP was associated with higher UV exposure; c the C allele of the MSTN SNP was associated with better performance at shorter racing distance; d the T allele of the MC1R SNP was associated with the chestnut coat colour when homozygous; and e the C allele of the NTM SNP was associated with more race starts in Thoroughbred racehorses

Discussion

High-resolution population structure of FM horses

Genetic forces (e.g. admixture and founder signatures) that have shaped the complex population structures of horses are currently inferred separately from each other using either parametric (e.g. Admixture) or non-parametric approaches (e.g. PCA). Recently, network-based cluster approaches, like NetView, are regaining favour for uncovering fine-scale population structures of horses [15, 17, 46]. The NetView approach can be used to visualize the relatedness of individuals without the need for any modelling of the genotype data [38]. Furthermore, it has been demonstrated that including individual admixture levels, genomic inbreeding, and genetic contributions in the network visualization describes a three-step procedure, which provides a high-resolution population structure without requiring prior ancestry information [35, 46].

Here, we applied the aforementioned three-step procedure to ascertain fine-scale population structure of FM horses. We showed that population substructures, as detected by model-based clustering, can be either associated with a different genetic origin or the progeny of most influential sires. For such complex population structures, we illustrated that identification of key contributors in combination with network visualization can shed light on the interpretation of admixture results, especially in situations where the cross-validation error estimation fails to provide an optimal number of clusters.

With regards to the formation of the FM breed, the high-resolution population structure analysis revealed that admixed FM horses showed unexpected high TB admixture levels, especially horses that descended from the Purebred Arabian stallion Doktryner. To clarify the origin of the DON line, further research is needed, including sequence analyses of the Y chromosome of male descendants [47,48,49,50]. Based on the identification of key contributors, we were also able to show that admixed FM horses (including one F1 sire from the N line and its progeny) had a strong impact on the current FM breeding population. Concerning the maintenance of genetic diversity in the breed, we further noticed that some FM individuals had lower FROH and fewer ROH segments (essentially the F1 and F2 crosses from the last introgressions), thereby decreasing the mean FROH. However, the frequency of very long ROH segments was high in both breeds (longer than 8 MB for more than 25% of segments). Both OF and FM appeared to have been genetically bottlenecked. For FM, it was mainly due to the intensive selection of stallions and breeding of rare coat colours (e.g. grey). In OF, the limitations imposed by the pure-breeding efforts has considerably restricted the gene pool. To prevent these inbreeding activities, we suggest to integrate OF horses in the FM breeding population and to encourage breeders to target traits other than rare coat colour, as selection for rare coat colour has been proven to be detrimental to maintaining genetic diversity [15].

Identification of ROH patterns and FROH in FM horses and other breeds

In our study, FROH for OF and FM was nearly twice as high than FPED for the majority of the horses, and higher than expected based on previous studies on FM using 50 K data, in which the mean FPED was slightly (0.2%) higher than FROH [51]. However, the medium correlation between FROH and FPED (r < 0.60) was similar to that in other livestock species such as goats [52], sheep [9], and cattle [53, 54]. In general, the mean FROH of a population tends to be overestimated compared to FPED: higher FROH than FPED have previously been reported in the Mangalarga Marchador [55], Italian Heavy Draught Horse [21], the Noriker, and to a much lesser extent in Haflinger horses [46]. Considering the depth of the pedigree data in the FM (> 20 generations), we suspect that both pedigree errors and selection of ROH settings affect the level of agreement between FROH and FPED. Currently, there are no standards for ROH settings in livestock [56] or horses. For the present study, we chose the same settings as previously applied in the study of other European horse breeds [14, 15, 17, 23]. We allowed for one heterozygous SNP per segment but phased the data to remove missing calls. However, as the genotype data is derived from SNP arrays, a further bias is introduced as the applied overlapping windows approach assumes that SNPs that are not included in the array are non-variants, which could cause an overestimation of ROH [56]. This will have little impact on the comparison of FROH between breeds in our study, but might affect the observed agreement (R2) between FROH and FPED of FM horses. However, by including information on parentage, the R2 between FROH and FPED increased considerably, as previously demonstrated [51]. Therefore, including parentage information in the estimation of FROH currently likely provides the most accurate information on inbreeding for conservation purposes.

With regards to other breeds, the TB had the most ROH segments, but they were also the shortest, indicative of historical inbreeding and consistent with an old established breed following a closed studbook [16, 57,58,59,60]. However, the cause for the higher FROH of TB compared to all other breeds may lie in the geographic origin of the samples, considering the strong bottleneck in the founder population of Australian TB [57]. The SA had the second shortest mean ROH segment lengths, and the third most ROH segments, indicative of a slightly less inbred population than TB. The AR had significantly more ROH segments, but they were not significantly longer than those of the SA population, although the AR had the highest proportion of very long ROH segments (> 10 Mb). The WB had ROH that were consistent with a breeding population with an open studbook, i.e. generally low inbreeding with some individual exceptions.

Runs of homozygosity islands and gene functions

We found no ROH islands for WB, in contrast to previous studies [13, 19, 21], which may be because we regrouped two European WB subpopulations (Swiss and Swedish) and by the selection of ROH settings that differed from the previous studies. The TB had many ROHs and these contained over 300 annotated genes (see Additional file 4: Table S3), including genes that were previously associated with memory formation (FSTL4), bone formation, lung function and development (WNT2B), muscle cell proliferation and regulation of myogenesis (KCNRG, MSTN), and hip conformation (IRX2) [61]. The ROH islands contained several interleukin genes (IL4, IL5, IL13), which, together with their associated receptors, have been implicated in inflammatory processes such as recurrent airway obstruction and summer-pasture associated obstructive pulmonary disease (IL4, IL13 [62,63,64]) or insect bite hypersensitivity (IL5 [65,66,67]). There were also three Insulin-like growth factor binding protein genes (IGFBP1, IGFBP3, and IGFBP6), which are involved in somatic cell growth and muscle cell distribution, and were previously identified in breed-specific ROH islands including German and Swedish WB horses [13, 19] as well as French Trotters [14].

The ROH islands for the two Arabian breeds (AR and SA) were previously described by Grilz-Seger et al. [14] and will therefore be only briefly discussed, although there were slight differences in the length and position of these ROH islands compared to this previous study due to the remapping of the SNP positions using the new EquCab3.0 reference and to the specific subset in the AR population investigated here (only Doktryner descendants) (see Additional file 2: Table S1 and Additional file 3: Table S2). The presence of a large cluster of Keratin genes, which are associated with consistency of skin, hair, and structures such as the hoof capsule, remains notable in both Arabian breeds. The GO term epithelial cell differentiation (GO:0030855) was in fact enriched in AR. A recent study also identified ZFMP1, a gene within an AR-specific ROH island, as a candidate gene for hoof health in Mongolian horses [68]. Overall, the presence of multiple genes involved in hoof structure suggests that the identified ROH islands may contribute to the overall hoof health of AR. Some genes found in AR- or TB-specific ROH islands were associated with performance: CBLB is considered a negative regulator of the insulin response to exercise and was identified in a transcriptomic study of Arabian horses in training [69]. The gene PPARGC1B was present in islands for AR and TB and was previously associated with endurance in horses [70] and humans [71]. Unsurprisingly, SA and AR shared several ROH islands, considering that SA was classified as a variety of Arabian until 1978 [72].

The well-characterised MSTN gene was present in the ROH islands for TB. The high frequency of the C-allele for MSTN in TB is consistent with the popularity of shorter races that are common in Australasia, while the oriental breeds (SA and AR especially) are mostly used in endurance racing and had a high frequency of the T-allele. This SNP was also associated with the distribution of muscle fibre types (C alleles for type IIA and IIX, T alleles for type I) in TB, but not in Belgian Draught horses, due to differences in linkage disequilibrium [73]. Generally, TB have higher proportions of fast twitch type IIA muscle fibres that are advantageous for anaerobic activity over a distance of 1000 m, while AR have high amounts of type I muscles fibres, which are used for aerobic endurance work [74]. It was expected that FM would have more T alleles and therefore more type I muscle fibres, which are solicited in draught training [75], similar to other draught breed such as the Belgian draught horse [73]. However, this hypothesis needs to be confirmed by muscle biopsies, as the exact interplay between the MSTN gene and muscle fibre type distribution has not been completely elucidated [73], and muscle fibre type distribution has not yet been studied in FM.

The FM had only few ROH islands compared to AR, SA, and TB (see Additional file 5: Table S4). The long ROH island for OF and FM on ECA11 (2.32 Mb), which was not shared with the other breeds, contained two genes that were previously identified in ROH islands in Noriker horses (LASP1) and in Noriker, Posavina, Lipizzaner, and Gidran horses (HOXB cluster) [14], as well as in German Warmblood horses [19]. Homeobox genes such as the HOXB complex are critical for development of the spine and appendages during embryological development [76]. It is therefore not surprising that the GO term embryonic skeletal system morphogenesis (GO:0048704) was enriched in OF and FM, while embryonic skeletal system development (GO:0048706) and skeletal system morphogenesis (GO:0048705) were additionally enriched in FM. Genes from the HOXB cluster were already under selection pressure in Europe during the Sassanid Persian expansion between the 7th to 9th century [77]. The FM, Noriker, Posavina, and Lipizzaner horses tend to favour a short back in relation to limb length, which might be the reason that  the HOXB gene cluster became relatively fixed in these breeds. However, without a definitive phenotype-genotype association, this remains speculative. The NOG gene, found in an ROH island for SA and FM but not in OF, is also associated with the length of phalanges [78, 79], which would suggest that breeders are selecting against shorter legs in FM compared to OF. The LASP1 gene has previously been identified as one of the major genes affecting size of the horse, particularly wither’s height [43]. However, based on the allele frequencies of the SNP that was previously identified by Makvandi-Nejad et al. [43], we concur with Grilz-Seger et al. [15] that the LASP1 gene is probably involved in body size (heavy or light type) rather than height itself. The FM had the highest frequency of “big” alleles, despite height at withers being limited to 160 cm in the breeding standards. Meanwhile, WB horses, who at least in Switzerland have to exceed 160 cm at the withers to be admitted to the studbook, have a low frequency of “big” alleles. However, SA had the lowest frequency of “big” alleles, although the SA should have a larger frame than the AR, and the “big” allele is even more frequent in FM than in OF [72]. This discrepancy is likely explained by other variants that did not appear in the ROH islands, such as ZFAT [15].

In an OF-specific ROH island (see Additional file 6 Table S5), we identified the gene Serotonin receptor 1B (HTR1B or 5-HT1B) on ECA10. Serotonin is an important neurotransmitter associated with coping styles in humans and animals [80]. HTR1B-knockout mice have been shown to be more reactive and less anxious than the wild-type mice [81] and more reactive to rewards [82]. Knockout rats are less likely to be startled but will react more forcefully under prolonged stress, indicative of a reactive coping style [83]. In horses, more reactive horses show no obvious reactions to aversive stimuli, which corresponds to the animals being known as “cold-blooded” [84, 85]. It is possible that the HTR1B gene has been under selection in OF as a working draught horse, for which strong flight-or-fight responses are undesirable. As this ROH island was not shared with FM, this would also imply that introgressions and selection towards a lighter type of horse influenced the genetic makeup of the breed as well as behaviour, towards a more proactive horse. However, HTR1B-knockout mice were also associated with an increase in aggressive behaviour [86], so it is not clear how the homozygous haplotype of HTR1B in the OF affects their behaviour. Considering the complexity of behaviour, more studies are needed to better interpret these results.

We identified another ROH island that has been associated with behaviour: SA, AR, and TB shared an ROH island containing the genes NTM and OPCLM, which belong to the IgLON protein family that regulate development of neuronal projections [87]. While downregulation of OPCLM has been associated with schizophrenia in Han Chinese [88], certain NTM polymorphisms in humans were associated with IQ [89] or aggressiveness in childhood ADHD [90]. NTM knockout mice were less fearful, moved more and reacted less to averse stimuli [91, 92], which might be favourable traits in sport horses, while OPCML knockout mice were less curious, had less spatial recognition and recognition memory in general, and showed less interest in their conspecifics [88]. In horses, the NTM gene has been one of the key genes involved in domestication [93] and has been found in selection signatures for racehorses such as TB [94] and Purebred Arabian horses [14]. Furthermore, this region on ECA7 has also been divergently selected between draught and light type breeds [95]. To better understand the direction of selection in NTM, we selected a SNP on the 670 K array that was previously described by McGivney et al. [45] to be associated with the number of race starts in Thoroughbreds. It was suggested that the favourable allele was indicative of a calmer temperament in the horse, with these horses racing earlier in life and therefore more often [45]. In our study, the favourable C allele was overrepresented in OF and FM, and nearly absent in all other breeds, even TB. This concurs with the breeding goal for a horse with good character. However, it could also be shown (see Additional file 12 Figure S2) that the unfavourable allele is associated with higher APED, suggesting that selection towards a lighter sport horse may affect the calm character the FM breed is known for.

Conclusions

The FM is the last native Swiss breed and is clearly distinguishable from other historically introgressed breeds, whereas admixed FM horses had a strong impact on the current breeding population. The FM was less inbred than other purebred horse breeds, but suffered bottlenecks due to intensive selection of stallions, restrictive mating choices based on arbitrary definitions of pure breeding, and selection for rare coat colours. Genes embedded in ROH islands had known functions involved in conformation and behaviour, two traits highly valued by breeders. The unexpected high TB admixture proportions in FM need further investigation.

Availability of data and materials

Anonymised data are available for academic purposes on signing a material transfer agreement.

References

  1. Poncet P-A, Wermeille V. Le cheval des Franches-Montagnes à travers l’histoire. Porrentruy: Société jurassienne d'émulation; 2009.

  2. Poncet P-A, Pfister W, Muntwyler J, Glowatzki-Mullis ML, Gaillard C. Analysis of pedigree and conformation data to explain genetic variability of the horse breed Franches-Montagnes. J Anim Breed Genet. 2006;123:114–21.

    Article  CAS  PubMed  Google Scholar 

  3. Hasler H, Flury C, Menet S, Haase B, Leeb T, Simianer H, et al. Genetic diversity in an indigenous horse breed–implications for mating strategies and the control of future inbreeding. J Anim Breed Genet. 2011;128:394–406.

    Article  CAS  PubMed  Google Scholar 

  4. Gmel AI, Burren A, Neuditschko M. Estimates of genetic parameters for shape space data in Franches-Montagnes horses. Animals. 2022;12:2186.

    Article  PubMed  PubMed Central  Google Scholar 

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

  6. McQuillan R, Leutenegger A-L, 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 

  7. Purfield DC, Berry DP, McParland S, Bradley DG. Runs of homozygosity and population history in cattle. BMC Genet. 2012;13:70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Zhang Q, Guldbrandtsen B, Bosse M, Lund MS, Sahana G. Runs of homozygosity and distribution of functional variants in the cattle genome. BMC Genomics. 2015;16:542.

    Article  PubMed  PubMed Central  Google Scholar 

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

  10. Signer-Hasler H, Burren A, Ammann P, Drögemüller C, Flury C. Runs of homozygosity and signatures of selection: A comparison among eight local Swiss sheep breeds. Anim Genet. 2019;50:512–25.

    Article  CAS  PubMed  Google Scholar 

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

  12. Gmel AI, Guichard M, Dainat B, Williams GR, Eynard S, Vignal A, et al. Identification of runs of homozygosity in Western honey bees (Apis mellifera) using whole-genome sequencing data. Ecol Evol. 2023;13: e9723.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Ablondi M, Viklund Å, Lindgren G, Eriksson S, Mikko S. Signatures of selection in the genome of Swedish warmblood horses selected for sport performance. BMC Genomics. 2019;20:717.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Grilz-Seger G, Neuditschko M, Ricard A, Velie B, Lindgren G, Mesarič M, et al. Genome-wide homozygosity patterns and evidence for selection in a set of European and near eastern horse breeds. Genes. 2019;10:491.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Grilz-Seger G, Druml T, Neuditschko M, Mesarič M, Cotman M, Brem G. Analysis of ROH patterns in the Noriker horse breed reveals signatures of selection for coat color and body size. Anim Genet. 2019;50:334–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Fawcett JA, Sato F, Sakamoto T, Iwasaki WM, Tozaki T, Innan H. Genome-wide SNP analysis of Japanese Thoroughbred racehorses. PLoS ONE. 2019;14: e0218407.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Grilz-Seger G, Druml T, Neuditschko M, Dobretsberger M, Horna M, Brem G. High-resolution population structure and runs of homozygosity reveal the genetic architecture of complex traits in the Lipizzan horse. BMC Genomics. 2019;20:174.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Metzger J, Karwath M, Tonda R, Beltran S, Águeda L, Gut M, et al. Runs of homozygosity reveal signatures of positive selection for reproduction traits in breed and non-breed horses. BMC Genomics. 2015;16:764.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Nolte W, Thaller G, Kuehn C. Selection signatures in four German warmblood horse breeds: tracing breeding history in the modern sport horse. PLoS ONE. 2019;14: e0215913.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Chen C, Zhu B, Tang X, Chen B, Liu M, Gao N, et al. Genome-wide assessment of runs of homozygosity by whole-genome sequencing in diverse horse breeds worldwide. Genes. 2023;14:1211.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Mancin E, Ablondi M, Mantovani R, Pigozzi G, Sabbioni A, Sartori C. Genetic variability in the Italian heavy draught horse from pedigree data and genomic information. Animals. 2020;10:1310.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Ablondi M, Dadousis C, Vasini M, Eriksson S, Mikko S, Sabbioni A. Genetic diversity and signatures of selection in a native Italian horse breed based on SNP data. Animals. 2020;10:1005.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Grilz-Seger G, Mesarič M, Cotman M, Neuditschko M, Druml T, Brem G. Runs of homozygosity and population history of three horse breeds with small population size. J Equine Vet Sci. 2018;71:27–34.

    Article  Google Scholar 

  24. Szmatoła T, Gurgul A, Jasielczuk I, Oclon E, Ropka-Molik K, Stefaniuk-Szmukier M, et al. Assessment and distribution of runs of homozygosity in horse breeds representing different utility types. Animals. 2022;12:3293.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Criscione A, Mastrangelo S, D’Alessandro E, Tumino S, Di Gerlando R, Zumbo A, et al. Genome-wide survey on three local horse populations with a focus on runs of homozygosity pattern. J Anim Breed Genet. 2022;139:540–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Drögemüller M, Jagannathan V, Welle MM, Graubner C, Straub R, Gerber V, et al. Congenital hepatic fibrosis in the Franches-Montagnes horse is associated with the polycystic kidney and hepatic disease 1 (PKHD1) gene. PLoS ONE. 2014;9: e110125.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Schaefer RJ, Schubert M, Bailey E, Bannasch DL, Barrey E, Bar-Gal GK, et al. Developing a 670k genotyping array to tag~ 2M SNPs across 24 horse breeds. BMC Genomics. 2017;18:565.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Kalbfleisch TS, Rice ES, DePriest MS, Walenz BP, Hestand MS, Vermeesch JR, et al. Improved reference genome for the domestic horse increases assembly contiguity and composition. Commun Biol. 2018;1:197.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Browning BL, Zhou Y, Browning SR. A one-penny imputed genome from next-generation reference panels. Am J Hum Genet. 2018;103:338–48.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience. 2015;4:7.

    Article  PubMed  PubMed Central  Google Scholar 

  31. R Core Team. R: A language and environment for statistical computing. Version 3.6.1. R Foundation for Statistical Computing; 2013.

  32. Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23:254–67.

    Article  CAS  PubMed  Google Scholar 

  33. Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19:1655–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8.

    Article  Google Scholar 

  35. Neuditschko M, Raadsma HW, Khatkar MS, Jonas E, Steinig EJ, Flury C, et al. Identification of key contributors in complex population structures. PLoS ONE. 2017;12: e0177638.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Dinno A, Dinno MA. Package ‘paran’. Version 1.5.2: CRAN (The Comprehensive R Archive Network); 2018.https://cran.r-project.org/web/packages/paran/index.html. Accessed 24 Apr 2024.

  37. Glorfeld LW. An improvement on Horn’s parallel analysis methodology for selecting the correct number of factors to retain. Educ Psychol Meas. 1995;55:377–93.

    Article  Google Scholar 

  38. Neuditschko M, Khatkar MS, Raadsma HW. NetView: a high-definition network-visualization approach to detect fine-scale population structures from genome-wide patterns of variation. PLoS ONE. 2012;7: e48375.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Steinig EJ, Neuditschko M, Khatkar MS, Raadsma HW, Zenger KR. netview p: a network visualization tool to unravel complex population structure using genome-wide SNPs. Mol Ecol Resour. 2016;16:216–27.

    Article  CAS  PubMed  Google Scholar 

  40. Graves S, Piepho H-P, Selzer ML. Package ‘multcompView’. Version 0.1-9: CRAN (The Comprehensive R Archive Network). 2015. https://cran.stat.unipd.it/web/packages/multcompView/multcompView.pdf. Accessed 24 Apr 2024.

  41. Coster A. Pedigree: package to deal with pedigree data. Version 1.4.2: CRAN (The Comprehensive R Archive Network). 2022. https://cran.r-project.org/web/packages/pedigree/index.html Accessed 11 June 2024.

  42. Biscarini F, Cozzi P, Gaspa G, Marras G. detectRUNS: an R package to detect runs of homozygosity and heterozygosity in diploid genomes. Version 0.9.6: CRAN (The Comprehensive R Archive Network). 2019. https://cran.r-project.org/web/packages/detectRUNS/index.html. Accessed 24 Apr 2024.

  43. Makvandi-Nejad S, Hoffman GE, Allen JJ, Chu E, Gu E, Chandler AM, et al. Four loci explain 83% of size variation in the horse. PLoS ONE. 2012;7: e39929.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Hill EW, Gu J, Eivers SS, Fonseca RG, McGivney BA, Govindarajan P, et al. A sequence polymorphism in MSTN predicts sprinting ability and racing stamina in thoroughbred horses. PLoS ONE. 2010;5: e8645.

    Article  PubMed  PubMed Central  Google Scholar 

  45. McGivney B, Hernandez B, Katz L, MacHugh D, McGovern S, Parnell A, et al. A genomic prediction model for racecourse starts in the Thoroughbred horse. Anim Genet. 2019;50:347–57.

    Article  CAS  PubMed  Google Scholar 

  46. Druml T, Neuditschko M, Grilz-Seger G, Horna M, Ricard A, Mesarič M, et al. Population networks associated with runs of homozygosity reveal new insights into the breeding history of the Haflinger horse. J Hered. 2018;109:384–92.

    Article  PubMed  Google Scholar 

  47. Bozlak E, Radovic L, Remer V, Rigler D, Allen L, Brem G, et al. Refining the evolutionary tree of the horse Y chromosome. Sci Rep. 2023;13:8954.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Felkel S, Vogl C, Rigler D, Dobretsberger V, Chowdhary BP, Distl O, et al. The horse Y chromosome as an informative marker for tracing sire lines. Sci Rep. 2019;9:6095.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Remer V, Bozlak E, Felkel S, Radovic L, Rigler D, Grilz-Seger G, et al. Y-Chromosomal insights into breeding history and sire line genealogies of Arabian horses. Genes. 2022;13:229.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Wallner B, Palmieri N, Vogl C, Rigler D, Bozlak E, Druml T, et al. Y chromosome uncovers the recent oriental origin of modern stallions. Curr Biol. 2017;27:2029-2035.e5.

    Article  CAS  PubMed  Google Scholar 

  51. Burren A, Signer-Hasler H, Neuditschko M, Leeb T, Rieder S, Flury C. Inzucht beim Freiberger Pferd. Agroscope Sci. 2016;32:38–9.

    Google Scholar 

  52. Burren A, Neuditschko M, Signer-Hasler H, Frischknecht M, Reber I, Menzi F, et al. Genetic diversity analyses reveal first insights into breed-specific selection signatures within Swiss goat breeds. Anim Genet. 2016;47:727–39.

    Article  CAS  PubMed  Google Scholar 

  53. Solé M, Gori A-S, Faux P, Bertrand A, Farnir F, Gautier M, et al. Age-based partitioning of individual genomic inbreeding levels in Belgian Blue cattle. Genet Sel Evol. 2017;49:92.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Ferenčaković M, Hamzić E, Gredler B, Solberg T, Klemetsdal G, Curik I, et al. Estimates of autozygosity derived from runs of homozygosity: empirical evidence from selected cattle populations. J Anim Breed Genet. 2013;130:286–93.

    Article  PubMed  Google Scholar 

  55. de Santos WB, Pimenta Schettini G, Fonseca MG, Pereira GL, Loyola Chardulo LA, Rodrigues Machado Neto O, 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 

  56. Lavanchy E, Goudet J. Effect of reduced genomic representation on using runs of homozygosity for inbreeding characterization. Mol Ecol Resour. 2023;23:787–802.

    Article  CAS  PubMed  Google Scholar 

  57. Todd ET, Ho SY, Thomson PC, Ang RA, Velie BD, Hamilton NA. Founder-specific inbreeding depression affects racing performance in Thoroughbred horses. Sci Rep. 2018;8:6167.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Orlando L, Librado P. Origin and evolution of deleterious mutations in horses. Genes. 2019;10:649.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Hill EW, McGivney BA, MacHugh DE. Inbreeding depression and durability in the North American Thoroughbred horse. Anim Genet. 2023;54:408–11.

    Article  PubMed  Google Scholar 

  60. Hill EW, Stoffel MA, McGivney BA, MacHugh DE, Pemberton JM. Inbreeding depression and the probability of racing in the Thoroughbred horse. Proc Biol Sci. 2022;289:20220487.

    CAS  PubMed  PubMed Central  Google Scholar 

  61. McGivney BA, Han H, Corduff LR, Katz LM, Tozaki T, MacHugh DE, et al. Genomic inbreeding trends, influential sire lines and selection in the global Thoroughbred horse population. Sci Rep. 2020;10:466.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Beadle R, Horohov D, Gaunt S. Interleukin-4 and interferon-gamma gene expression in summer pasture-associated obstructive pulmonary disease affected horses. Equine Vet J. 2002;34:389–94.

    Article  CAS  PubMed  Google Scholar 

  63. Klukowska-Rötzler J, Swinburne J, Drögemüller C, Dolf G, Janda J, Leeb T, et al. The interleukin 4 receptor gene and its role in recurrent airway obstruction in Swiss Warmblood horses. Anim Genet. 2012;43:450–3.

    Article  PubMed  Google Scholar 

  64. Korn A, Miller D, Dong L, Buckles EL, Wagner B, Ainsworth DM. Differential gene expression profiles and selected cytokine protein analysis of mediastinal lymph nodes of horses with chronic recurrent airway obstruction (RAO) support an interleukin-17 immune response. PLoS ONE. 2015;10: e0142622.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Langreder N, Schäckermann D, Meier D, Becker M, Schubert M, Dübel S, et al. Development of an inhibiting antibody against equine interleukin 5 to treat insect bite hypersensitivity of horses. Sci Rep. 2023;13:4029.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Jonsdottir S, Fettelschoss V, Olomski F, Talker SC, Mirkovitch J, Rhiner T, et al. Safety profile of a virus-like particle-based vaccine targeting self-protein interleukin-5 in horses. Vaccines. 2020;8:213.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Pantelyushin S, Rhiner T, Jebbawi F, Sella F, Waldern N, Lam J, et al. Interleukin 5-dependent inflammatory eosinophil subtype involved in allergic insect bite hypersensitivity of horses. Allergy. 2023;78:3020–3.

    Article  CAS  PubMed  Google Scholar 

  68. Han H, Randhawa IA, MacHugh DE, McGivney BA, Katz LM, Dugarjaviin M, et al. Selection signatures for local and regional adaptation in Chinese Mongolian horse breeds reveal candidate genes for hoof health. BMC Genomics. 2023;24:35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Ropka-Molik K, Stefaniuk-Szmukier M, Żukowski K, Piórkowska K, Gurgul A, Bugno-Poniewierska M. Transcriptome profiling of Arabian horse blood during training regimens. BMC Genet. 2017;18:31.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Mach N, Plancade S, Pacholewska A, Lecardonnel J, Rivière J, Moroldo M, et al. Integrated mRNA and miRNA expression profiling in blood reveals candidate biomarkers associated with endurance exercise in the horse. Sci Rep. 2016;6:22932.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Akhmetov I, Popov D, Missina S, Vinogradova O, Rogozkin V. The analysis of PPARGC1B gene polymorphism in athletes. Ross Fiziol Zh Im IM Sechenova. 2009;95:1247–53.

    CAS  Google Scholar 

  72. Doliş M, Nagy P, Doliş L, Nistor C. Study on the evolution of the body height sizes in female young horses of the Shagya breed. Lucr Științ Univ Agron Ser Zooteh Med Vet. 2011;56:219–23.

    Google Scholar 

  73. Petersen JL, Valberg SJ, Mickelson JR, McCue ME. Haplotype diversity in the equine myostatin gene with focus on variants associated with race distance propensity and muscle fiber type proportions. Anim Genet. 2014;45:827–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. López-Rivero JL, Agüera E, Monterde JG, Rodríguez-Barbudo MV, Miró F. Comparative study of muscle fiber type composition in the middle gluteal muscle of Andalusian, Thoroughbred and Arabian horses. J Equine Vet Sci. 1989;9:337–40.

    Article  Google Scholar 

  75. Serrano AL, Rivero JLL. Myosin heavy chain profile of equine gluteus medius muscle following prolonged draught-exercise training and detraining. J Muscle Res Cell Mot. 2000;21:235–45.

    Article  CAS  Google Scholar 

  76. Wellik DM. Hox patterning of the vertebrate axial skeleton. Dev Dyn. 2007;236:2454–63.

    Article  CAS  PubMed  Google Scholar 

  77. Fages A, Hanghøj K, Khan N, Gaunitz C, Seguin-Orlando A, Leonardi M, et al. Tracking five millennia of horse management with extensive ancient genome time series. Cell. 2019;177:1419-35.e31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Khan S, Mudassir M, Khan N, Marwat A. Brachdactyly instigated as a result of mutation in GDF5 and NOG genes in Pakistani population. Pak J Med Sci. 2018;34:82–7.

    Article  PubMed  PubMed Central  Google Scholar 

  79. Lehmann K, Seemann P, Silan F, Goecke T, Irgang S, Kjaer K, et al. A new subtype of brachydactyly type B caused by point mutations in the bone morphogenetic protein antagonist NOGGIN. Am J Hum Genet. 2007;81:388–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Puglisi-Allegra S, Andolina D. Serotonin and stress coping. Behav Brain Res. 2015;277:58–67.

    Article  CAS  PubMed  Google Scholar 

  81. Zhuang X, Gross C, Santarelli L, Compan V, Trillat A-C, Hen R. Altered emotional states in knockout mice lacking 5-HT1A or 5-HT1B receptors. Neuropsychopharmacology. 1999;21:52S-60S.

    Article  CAS  PubMed  Google Scholar 

  82. Crabbe JC, Phillips TJ, Feller DJ, Hen R, Wenger CD, Lessov CN, et al. Elevated alcohol consumption in null mutant mice lacking 5–HT1B serotonin receptors. Nat Genet. 1996;14:98–101.

    Article  CAS  PubMed  Google Scholar 

  83. Clark MS, Vincow ES, Sexton TJ, Neumaier JF. Increased expression of 5-HT1B receptor in dorsal raphe nucleus decreases fear-potentiated startle in a stress dependent manner. Brain Res. 2004;1007:86–97.

    Article  CAS  PubMed  Google Scholar 

  84. Seaman S, Davidson H, Waran N. How reliable is temperament assessment in the domestic horse (Equus caballus)? Appl Anim Behav Sci. 2002;78:175–91.

    Article  Google Scholar 

  85. Briefer Freymond S, Bardou D, Beuret S, Bachmann I, Zuberbühler K, Briefer EF. Elevated sensitivity to tactile stimuli in stereotypic horses. Front Vet Sci. 2019;6:162.

    Article  PubMed  PubMed Central  Google Scholar 

  86. Saudou F, Amara DA, Dierich A, LeMeur M, Ramboz S, Segu L, et al. Enhanced aggressive behavior in mice lacking 5-HT1B receptor. Science. 1994;265:1875–8.

    Article  CAS  PubMed  Google Scholar 

  87. Hashimoto T, Maekawa S, Miyata S. IgLON cell adhesion molecules regulate synaptogenesis in hippocampal neurons. Cell Biochem Funct. 2009;27:496–8.

    Article  CAS  PubMed  Google Scholar 

  88. Zhang Z, Ye M, Li Q, You Y, Yu H, Ma Y, et al. The schizophrenia susceptibility gene OPCML regulates spine maturation and cognitive behaviors through Eph-Cofilin signaling. Cell Rep. 2019;29:49-61.e47.

    Article  CAS  PubMed  Google Scholar 

  89. Pan Y, Wang K-S, Aragam N. NTM and NR3C2 polymorphisms influencing intelligence: family-based association studies. Prog Neuro-Psychopharmacol Biol Psychiatry. 2011;35:154–60.

    Article  CAS  Google Scholar 

  90. Brevik EJ, van Donkelaar MM, Weber H, Sánchez-Mora C, Jacob C, Rivero O, et al. Genome-wide analyses of aggressiveness in attention-deficit hyperactivity disorder. Am J Med Genet B Neuropsychiatr Genet. 2016;171:733–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. Singh K, Lilleväli K, Gilbert SF, Bregin A, Jayaram M, Rahi M, et al. The combined impact of IgLON family proteins Lsamp and Neurotrimin on developing neurons and behavioral profiles in mouse. Brain Res Bull. 2018;140:5–18.

    Article  CAS  PubMed  Google Scholar 

  92. Mazitov T, Bregin A, Philips M-A, Innos J, Vasar E. Deficit in emotional learning in neurotrimin knockout mice. Behav Brain Res. 2017;317:311–8.

    Article  CAS  PubMed  Google Scholar 

  93. Schubert M, Jónsson H, Chang D, Der Sarkissian C, Ermini L, Ginolhac A, et al. Prehistoric genomes reveal the genetic foundation and cost of horse domestication. Proc Natl Acad Sci USA. 2014;111:E5661–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Han H, McGivney BA, Allen L, Bai D, Corduff LR, Davaakhuu G, et al. Common protein-coding variants influence the racing phenotype in galloping racehorse breeds. Commun Biol. 2022;5:1320.

    Article  PubMed  PubMed Central  Google Scholar 

  95. Gurgul A, Jasielczuk I, Semik-Gurgul E, Pawlina-Tyszko K, Stefaniuk-Szmukier M, Szmatoła T, et al. A genome-wide scan for diversifying selection signatures in selected horse breeds. PLoS ONE. 2019;14: e0210751.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We gratefully acknowledge the support from the Swiss Franches-Montagnes breeding association for the provided pedigree information, and the Fondation Sur-La-Croix for funding part of the genotyping for the Franches-Montagnes horses.

Funding

This study was funded by the Swiss Federal Office for Agriculture (FOAG) under contract number 625000469 and Fondation Sur-La-Croix (internal contract number 6510263).

Author information

Authors and Affiliations

Authors

Contributions

AIG: conceptualisation, methodology, data collection, formal analysis, visualisation, writing—original draft. SM, AR, BV, VG: data collection and curation, writing—review and editing. MN: conceptualisation, methodology, software, resources, supervision, data curation, formal analysis, visualisation, writing—original draft. All authors have read and approved the manuscript.

Corresponding author

Correspondence to Markus Neuditschko.

Ethics declarations

Ethics approval and consent to participate

This study was conducted on data sampled under permits VD3527b, VD2976.1 and VD2227.2, each approved by the cantonal veterinary office of Vaud, Switzerland.

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

12711_2024_922_MOESM1_ESM.pdf

Additional file 1: Figure S1. Identification of the optimal number of clusters (K). Illustration of the cross-validation error increasing K from 1 to 20.

12711_2024_922_MOESM2_ESM.pdf

Additional file 2: Table S1. Runs of homozygosity segments shared by more than 50% of Purebred Arabians. Table S1 presents the runs of homozygosity segments shared by more than 50% of Purebred Arabians including the length and position along the chromosomes. The annotated genes within the segments are also reported.

12711_2024_922_MOESM3_ESM.pdf

Additional file 3: Table S2. Runs of homozygosity segments shared by more than 50% of Shagya Arabians. Table S2 presents the runs of homozygosity segments shared by more than 50% of Shagya Arabians including the length and position along the chromosomes. The annotated genes within the segments are also reported.

12711_2024_922_MOESM4_ESM.pdf

Additional file 4: Table S3. Runs of homozygosity segments shared by more than 50% of Thoroughbreds. Table S3 presents the runs of homozygosity segments shared by more than 50% of Thoroughbreds including the length and position along the chromosomes. The annotated genes within the segments are also reported.

12711_2024_922_MOESM5_ESM.pdf

Additional file 5: Table S4. Runs of homozygosity segments shared by more than 50% of Old-type Franches-Montagnes. Table S4 presents the runs of homozygosity segments shared by more than 50% of Old-type Franches-Montagnes including the length and position along the chromosomes. The annotated genes within the segments are also reported.

12711_2024_922_MOESM6_ESM.pdf

Additional file 6: Table S5. Runs of homozygosity segments shared by more than 50% of Franches-Montagnes. Table S5 presents the runs of homozygosity segments shared by more than 50% of Franches-Montagnes including the length and position along the chromosomes. The annotated genes within the segments are also reported.

12711_2024_922_MOESM7_ESM.pdf

Additional file 7: Table S6. Gene ontology analysis for the Purebred Arabians. Table S6 presents the results from the gene ontology analysis for the genes present in the runs of homozygosity islands for the Purebred Arabians.

12711_2024_922_MOESM8_ESM.pdf

Additional file 8: Table S7. Gene ontology analysis for the Shagya Arabians. Table S7 presents the results from the gene ontology analysis for the genes present in the runs of homozygosity islands for the Shagya Arabians.

12711_2024_922_MOESM9_ESM.pdf

Additional file 9: Table S8. Gene ontology analysis for the Thoroughbred. Table S8 presents the results from the gene ontology analysis for the genes present in the runs of homozygosity islands for the Thoroughbred.

12711_2024_922_MOESM10_ESM.pdf

Additional file 10: Table S9. Gene ontology analysis for the old-type Franches-Montagnes. Table S9 presents the results from the gene ontology analysis for the genes present in the runs of homozygosity islands for the old-type Franches-Montagnes.

12711_2024_922_MOESM11_ESM.pdf

Additional file 11: Table S10. Gene ontology analysis for the Franches-Montagnes. Table 10 presents the results from the gene ontology analysis for the genes present in the runs of homozygosity islands for the Franches-Montagnes.

12711_2024_922_MOESM12_ESM.pdf

Additional file 12: Figure S2. Association between Percent of pedigree-based admixture proportion and NTM SNP genotype. The C allele of the NTM SNP was associated with more race starts in Thoroughbred racehorses. Boxplot representing the pedigree-based admixture based on the NTM SNP genotype.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Gmel, A.I., Mikko, S., Ricard, A. et al. Using high-density SNP data to unravel the origin of the Franches-Montagnes horse breed. Genet Sel Evol 56, 53 (2024). https://doi.org/10.1186/s12711-024-00922-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12711-024-00922-6