High-resolution population structure and runs of homozygosity reveal the genetic architecture of complex traits in the Lipizzan horse

Background The sample ascertainment bias due to complex population structures remains a major challenge in genome-wide investigations of complex traits. In this study we derived the high-resolution population structure and levels of autozygosity of 377 Lipizzan horses originating from five different European stud farms utilizing the SNP genotype information of the high density 700 k Affymetrix Axiom™ Equine genotyping array. Scanning the genome for overlapping runs of homozygosity (ROH) shared by more than 50% of horses, we identified homozygous regions (ROH islands) in order to investigate the gene content of those candidate regions by gene ontology and enrichment analyses. Results The high-resolution population network approach revealed well-defined substructures according to the origin of the horses (Austria, Slovakia, Croatia and Hungary). The highest mean genome coverage of ROH (SROH) was identified in the Austrian (SROH = 342.9), followed by Croatian (SROH = 214.7), Slovakian (SROH = 205.1) and Hungarian (SROH = 171.5) subpopulations. ROH island analysis revealed five common islands on ECA11 and ECA14, hereby confirming a closer genetic relationship between the Hungarian and Croatian as well as between the Austrian and Slovakian samples. Private islands were detected for the Hungarian and the Austrian Lipizzan subpopulations. All subpopulations shared a homozygous region on ECA11, nearly identical in position and length containing among other genes the homeobox-B cluster, which was also significantly (p < 0.001) highlighted by enrichment analysis. Gene ontology terms were mostly related to biological processes involved in embryonic morphogenesis and anterior/posterior specification. Around the STX17 gene (causative for greying), we identified a ROH island harbouring the genes NR4A3, STX17, ERP44 and INVS. Within further islands on ECA14, ECA16 and ECA20 we detected the genes SPRY4, NDFIP1, IMPDH2, HSP90AB1, whereas SPRY4 and HSP90AB1 are involved in melanoma metastasis and survival rate of melanoma patients in humans. Conclusions We demonstrated that the assessment of high-resolution population structures within one single breed supports the downstream genetic analyses (e.g. the identification of ROH islands). By means of ROH island analyses, we identified the genes SPRY4, NDFIP1, IMPDH2, HSP90AB1, which might play an important role for further studies on equine melanoma. Furthermore, our results highlighted the impact of the homeobox-A and B cluster involved in morphogenesis of Lipizzan horses. Electronic supplementary material The online version of this article (10.1186/s12864-019-5564-x) contains supplementary material, which is available to authorized users.


Background
The Lipizzan horse breed is globally one of the best-documented horse population, as pedigree records can be traced back to the known founder animals born in the early eighteenth century. The founder population, described in detail by Zechner et al. [1], Druml and Sölkner [2] and Druml et al. [3], comprises 456 animals, whereas the major part of horses originated from the former imperial stud farm of Lipica founded in the year 1580. Since the First World War the stud book of Lipizzan horses is closed and presently conservation breeding strategies are applied by eleven state stud farms located in nine European countries. National Lipizzan breeding herds have limited population sizes; therefore, maintenance of genetic diversity has been in the focus of state stud farms over decennia.
Numerous genetic approaches and methods have been applied to characterise the gene pool, the genetic diversity and the population structure of the Lipizzan breed. Comprehensive pedigree data (max. 31 generations, generation equivalent 19.1) were used to estimate inbreeding coefficients and effective population sizes by Zechner et al. [1], Sölkner and Druml [2] and Druml et al. [3]. Furthermore, microsatellite markers were employed to describe diversity measures, genetic distances and population structure by Achmann et al. [4]. Kavar et al. [5,6] investigated mtDNA maternal diversity; Kasarda et al. [7] estimated genetic relatedness between Old Kladruber, Slovenian and Slovakian Lipizzans, whilst Wallner et al. [8] highlighted the patrilinear structure conducting haplotype-based analyses of the Y-chromosome. Most of these scientific publications arose from a multilateral research project, which is described in the review by Dovc et al. [9]. Based upon this scientific project further research was conducted focusing on the inheritance of melanoma, vitiligo and greying [10][11][12][13].
In this study we used the SNP genotype information of the Affymetrix Axiom™ Equine genotyping array [14] to analyse the high-resolution population structure of Lipizzan horses originating from five different European stud farms, which differ in breeding history and breeding objectives. To ascertain the high-resolution population structure of the gene pool, we applied a recently described three-step procedure as presented in Druml et al. [15], which includes individual levels of admixture and genomic inbreeding (ROH) in the final population network visualization. After the assessment of the high-resolution population structure, we identified overlapping homozygous regions (ROH islands) within the entire population and respective subpopulations, and conducted a gene ontology and enrichment analysis of annotated genes located within ROH islands. We demonstrate that the combination of different approaches elaborate insights in population structure and underlying differences in levels of autozygosity and distribution of ROH islands.

Genetic diversity and admixture analysis
Principal Component analysis (PCA) as presented in Fig. 1a) revealed that the first three Principal Components (PCs) accounted for 35% of the total genetic variance. Based upon the visualization of the first three PCs the Austrian population is characterized by lower pairwise genetic distances and a distinct clustering. The other three subpopulations (Croatia, Hungary and Slovakia) appear to be highly interrelated simultaneously expressing a high level of genetic diversity, whilst a subset of the Austrian samples shares similarities with Croatian, Slovakian and Hungarian horses. F ST analysis recapitulated these findings, as nearest relationships were documented for Austrian, Croatian and Slovakian samples (pairwise F ST from 0.025 to 0.045), whereas the Hungarian subpopulation exhibited higher genetic differentiation (pairwise F ST from 0.068 to 0.092) (Additional file 1).
The first level (K = 2) of model-based clustering using the programme Admixture separated the Austrian population from the other Lipizzan samples (Fig. 1b). At the second level of clustering K = 3 the Hungarian sample were allocated in a distinct cluster. Further increasing K to 4 and 5 the Austrian subpopulation was further sub-structured, whilst at K = 6 the Slovakian sample formed a distinct cluster. At the additional levels of K = 7 and K = 8 the identified subpopulations were subsequently sub-structured. The visualization of the cross-validation (CV) error for each K, increasing K from 2 to 10, did not result in an optimal number of clusters (Additional file 2).

High-resolution network visualisation
Using a high-resolution network graph we combined individual levels of S ROH (see the results below) and individual levels of admixture (K = 8). The network analysis improved the assessment of the population structure, compared to PCA and Admixture, by clearly separating the horses according to their origin, simultaneously highlighting four population outliers (Fig. 2). The horses of the stud farm Piber characteristically were connected by short and thick links, representing higher co-ancestry and lower genetic distance. The network highlighted two sub-clusters with higher levels of admixture, which were linked to Slovakian and Hungarian Lipizzans (marked with arrows in Fig. 2). According to pedigree information these cross-link animals were foreign bred mares, which were integrated into the Piber breeding population and thus are characterized by a lower co-ancestry level (longer and thinner links) and smaller S ROH values (smaller node sizes).  The Croatian sample was separated from the Austrian sample, connected directly by one cross-link animal. The network analysis identified three sub-clusters within the Croatian sample. According to pedigree information, the cluster on the right involved horses from the stud farm Lipik, whilst the left cluster represented the core population from the stud farm Đakovo. The accumulation in the middle of the Croatian sample included horses from Đakovo, which showed high levels of admixture with horses from Lipik, as well as with Slovakian and Hungarian Lipizzans, respectively. The Hungarian sample was characterised by the highest genetic distances between the samples and by the lowest S ROH values, indicating genealogical separation between Hungarian and Austrian Lipizzans. The Slovakian sample showed genetic relationships in both directionsto the Austrian and to the Croatian samples. These relationships were illustrated by manifold cross-link animals, which had small S ROH values and higher degree of admixture connecting both clusters.

Runs of homozygosity
Including the entire sample of Lipizzan horses in the ROH analysis the mean number of ROH (N ROH ) comprised 202.07 segments at a mean genome length covered by ROH (S ROH ) of 297.01 Mb (+ − 119.85), resulting in a mean ROH length (L ROH ) of 1.42 Mb (Table 1) Inbreeding as defined by F ROH reached a mean value of 13.2% (+ − 5.3%) within the entire sample, ranging from 15.3% (+ − 4.7%) in the Piber population to 7.6% (+ − 5.0%) in the Hungarian sample (Table 1). To differentiate old and recent inbreeding we calculated F ROH based upon different ROH length classes ( Table 2). The values for F ROH < 2Mb were highest in the Austrian (6.7%) and lowest in the Hungarian sample (5.6%). To quantify recent inbreeding we calculated F ROH considering ROHs longer than 10 MB. Recent inbreeding was absent in the Slovakian sample and tended to zero in the Hungarian sample, whereas the highest value of 1.1% was observed for the Austrian sample.
ROH patterns as revealed by the plot N ROH versus S ROH indicate for the majority of the Piber population a strong right shift towards a higher S ROH and an increased variance of S ROH (Additional file 3). The samples of Slovakia and Croatia were distributed along the diagonal, whereas in the left corner of the plot a cluster of 11 horses with the lowest ROH parameters can be (See figure on previous page.) Fig. 2 High-resolution network visualisation of Lipizzan horses from the stud farms Piber (Austria), Topol'čianky (Slovakia), Đakovo and Lipik (Croatia) and Szilvasvárad (Hungary). a On the top network the horses (nodes) are colored by stud farm. On the lower network (b) nodes are containing individual cluster membership coefficients from Admixture analysis. Each node represents an individual and the size of node is in relation to S ROH , link lengths between nodes illustrate the genetic distance between individuals, thickness of links are proportional to genetic relationship. Two arrows pinpoint sub-clusters in the Austrian sample observed, indicating out-breeding. These horses were also clearly identified as cross-link animals between stud farms within the high-resolution population network (Fig. 2). Concerning the distribution of ROH segments of different length categories, the Austrian Lipizzan population is characterized by the lowest proportion of ROH smaller than 2 Mb (21.3%), which varied in the other samples between 31.1 and 32.3% (Table 2). At the same time the Austrian horses had the highest proportion of ROHs longer than 6 Mb (19.2%). In the other samples this ROH category was present at an amount from 1.4 to 3.2%. In the Austrian sample 7.4% of ROHs were longer than 10 Mb, whereas in the other samples this category was underrepresented at a percentage of < 0.6%.

ROH islands
Within the entire dataset, we identified three ROH islands on ECA11 and two islands on ECA14, which were shared by more than 50% of the horses ( Table 3). The overlapping homozygous region on ECA11:24.13-24.81 Mb was shared by 61.3% of animals and harboured the genes COPZ2, CBX1, SNX11, HOXB1, HOXB2, HOXB3, HOXB5, HOXB6, HOXB7, HOXB8, HOXB13, TTLL6. Figure 3a illustrates this island and the position of the homeobox-B cluster (HOXB) for each sample separately. These overlapping homozygous regions were nearly identical in position and length in all samples and the frequency varied between 59.2% (Slovakia) and 72.7% (Hungary).
The second ROH island on ECA11 (position 31,015.821 -31,942.748) was up to 927 kb long in Austrian and Croatian samples and was shared by 72.8 and 62.6% of the horses, respectively (Fig. 3b). The Slovakian Lipizzans had two overlapping homozygous segments in this region. Within the Hungarian sample, 47.0% of horses shared a much smaller island, which was directly located in the Musashi RNA binding protein 2 (MSI2).
We identified a ROH island on ECA25 around the STX17 (syntaxin-17) gene responsible for grey coat colour, shared by 46.2% of horses within the entire sample, containing the genes NR4A3, STX17, ERP44 and INVS. This ROH island was embedded in the centre of a 1.38 Mb long homozygous region at position ECA25:5.69-7.07 Mb present in 36.6% of Lipizzan horses. The frequency was below our expectations. From a previous study, it is known that the STX17 genotype frequencies for homozygous G/G horses reached up to 67.3% in Lipizzans [13]. Scanning this region for ROHs shorter than 500 kb, we applied a smaller boundary (80 kb window length, min. 20 SNPs), and detected a 399.9kb long ROH island at position ECA25:6394.110-6794.044 centred round the STX17 gene that was shared by 66.1% of the entire sample ( Fig. 4). Across the subpopulations the frequencies for this ROH   island was highest in the Austrian sample (71.2%) and lowest in the Croatian sample (51.1%).
Within the Slovakian sample we detected six ROH islands on ECA11 (3), ECA14, ECA16 and ECA22, whereas the islands on ECA11 and ECA16 overlapped with corresponding islands of the Austrian Lipizzans ( Table 5). The ROH islands on ECA22:47.90-48.48 Mb and ECA14:34.65-35.08 Mb were also identified in the Croatian sample.
Within the Croatian sample six ROH islands on ECA4, ECA7, ECA11 (2 islands), ECA14, ECA22, which varied in length between 129.1 kb and 1.14 Mb, were identified. The longest overlapping homozygous region on ECA14:34.05-35.18 Mb was shared by 72.6% of horses and contained nine annotated genes ( Table 6).
The island on ECA4 at position 58.11-58.62 Mb was also present in the Hungarian Lipizzans and contained among other annotated genes the homeobox-A cluster (HOXA). Furthermore, an island on ECA22 at position 47.90-48.62 Mb was exclusively detected within the Croatian and Slovakian sample.
Within the Hungarian Lipizzans six ROH islands on ECA4, ECA8, ECA11, ECA22, ECA23, ECA30 were identified, varying in length between 497.4 kb and 699.2 kb (Table 7). Four ROH islands on ECA8, ECA22, ECA23 and ECA30 were specific for the Hungarian sample. On ECA8 also a private island of the Austrian sample was documented (Fig. 5).

Gene ontology and enrichment analysis
Gene ontology and enrichment analysis highlighted within all Lipizzan samples the homeobox-B cluster, where the highest significance levels (p < 0.001) were reached for the terms GO:0048704~embryonic skeletal system morphogenesis and GO:0009952~anterior/posterior pattern specification (Additional file 4). Further the term GO:0043565~sequence-specific DNA binding related to molecular functions reached a significance level of p < 0.001 for the entire Lipizzan sample. Three out of 20 genes located in the identified ROH islands were listed in OMIM (Online Mendelian Inheritance in Man) database (DGKE -615,008~Hemolytic uremic syndrome; HOXB1-614744~Facial paresis; SPRY4-615266~Hypogonadotropic hypogonadism), but none of the identified genes were listed in OMIA (Online Mendelian Inheritance in Animals) database for horses.
In the Austrian Lipizzan population (Additional file 5) GO analysis confirmed the aforementioned terms in biological processes (GO:0048704 and GO:0009952 containing the HOXB cluster), and highlighted additionally the five terms GO:0021570~rhombomere 4 development, GO:0021612~facial nerve structural, GO:0006183~GTP biosynthetic process, GO:0006950~response to stress and GO:0071353~cellular response to interleukin-4, related to biological processes.
The Slovakian Lipizzans (Additional file 6) where characterized by three specific terms related to biological processes: GO:0048864~stem cell development, GO:0001525~angiogenesis, GO:0006368~transcription elongation from RNA polymerase II promoter.
GO analysis revealed 12 identical terms, related to biological processes for the Hungarian and Croatian samples and additionally two sample specific terms (GO:0006355~regulation of transcription, DNA-templated; GO:0001578m icrotubule bundle formation) for the Hungarian and three terms for the Croatian sample (GO:0007156~homophilic cell adhesion via plasma membrane adhesion molecules; GO:0050890~cognition; GO:0001759~organ induction)  (Additional files 7 and 8). All these 12 terms related to biological processes are based upon the HOXB-and HOX-A-clusters. The highest significance levels (p < 0.001) were found in concordance with the other samples for GO:0009952~anterior/posterior pattern specification and GO:0048704~embryonic skeletal system morphogenesis. Summarizing the results from gene ontology analysis of these two samples a high number of terms were highlighted for different aspects of embryonic morphogenesis (f.e. GO:0009953~dorsal/ventral pattern formation, GO:0009954~proximal/distal pattern formation, GO:0008584~male gonad development, GO:0060065~uterus development, GO:0021615~glossopharyngeal nerve morphogenesis).

Single SNP analysis
Two islands in the entire sample and six islands in the Austrian Lipizzans were remarkable short (8.

Discussion
The high-resolution population network (Admixture and Netview analysis) of the Lipizzan horse breed revealed that the single stud farms represent differentiated subpopulations, whereas the highest genetic distances were observed between the Austrian and the Hungarian/Croatian samples. The Slovakian Lipizzans clustered between the Austrian and the Croatian Lipizzans. These findings are in concordance with Achmann et al. [4], who derived a comparable structure for 561 Lipizzans from seven European countries. In the present study the Austrian and Slovakian horses were characterized by a lower genetic distances and higher pairwise relationship coefficients compared to the Croatian and Hungarian animals. This can be explained by the limited census of breeding animals in the Slovakian stud farm Topol'čianky, and the limited level of introgression of horses into the Austrian breeding herd. Due to the obligatory performance test of stallions in the Spanish riding school, in the Austrian stud farm Piber only foreign mares can be used to increase genetic diversity, whereas in Hungary, Slovakia and Croatia an exchange of breeding stallions is commonly applied. At the same time the Croatian and Hungarian breeding stock has been bigger than the Slovakian one. Generally genetic distances between horses, illustrated in the high-resolution network visualization, was relatively equally distributed within the breeding herds of Piber and Topol'čianky and no outstanding contribution of single animals (key-contributor) were detected, revealing the efforts of a conservation breeding program with narrow sex ratios and moderate selection intensity. Between Piber, Szilvasvárad and Topol'čianky cross-link animals were found, which showed markedly smaller mean genome-wide coverage of ROH and a higher degree of admixture. A comparable "outbreeding-effect" was detected in a previous study in the Haflinger breed, where outcrosses between horses from different countries/stud books were characterised by low S ROH and the lack of ROHs > 6 Mb [15].
All ROH parameters observed in our four subpopulations, reached highest values in the Austrian Lipizzans, whereas the lowest values were found in the Hungarian sample. The values of F ROH generally were in concordance with the findings of pedigree analyses from Zechner et al.  [1] and did not exceed 15%. Genomic inbreeding, as described by F ROH in the Lipizzan breed (F ROH = 0.13) was lower than expected from long-term stud farm breeding at small census and it is comparable to values found in the Austrian Haflinger (0.13), the Slovenian Haflinger (0.12) and the Bosnian Mountain Horse (0.14) [15,16]. Higher levels of F ROH were identified in Shagya Arabians (F ROH = 0.16) and Purebred Arabians (F ROH = 0.18), two breeds characterized by early closure of stud books [15].The length of ROHs is expected to correlate to ancient and recent inbreeding due to number of recombination events. Thus, recent inbreeding events result in longer ROHs as only a few identical by descent (IBD) segments can be broken down by repeated meiosis. According to Browning and Browning [17] and Thompson [18] ROHs of a length of 16.6, 10.0 and 5.0 Mb are assumed to origin from common ancestors back in the 3rd, 5th and 10th generations (6, 10 and 20 meioses respectively). Considering a mean generation interval of 10 years in Lipizzan horses the origin of ROHs in length class of < 2 Mb can be dated to the foundation period of this breed (foundation time 1580 AC). We observed F ROH values calculated based upon ROHs shorter than 2 Mb ranging from 5.6 to 7.1%. F ROH > 10Mb was absent in Slovakian Lipizzans and tended toward zero in Croatian and Hungarian horses, indicating minimal recent inbreeding within the past five generations [18,19]. In the Austrian population F ROH > 10Mb reached 1.1%. These results demonstrated that it is feasible to minimize the increase of inbreeding in populations limited in size and closed stud books. The ROH profile illustration for the Austrian Lipizzan revealed a strong right shift and a higher variance typical for consanguineous and bottlenecked populations  [20]. The Slovakian sample was characterized by a similar profile, but at a much lower extend, whereas the Croatian and Hungarian samples showed the effects of small census [20]. The identification of overlapping homozygous regions confirmed a closer genetic relationship between the Hungarian and Croatian samples as well as between the Austrian and Slovakian samples, which was also highlighted within the high-resolution population network. Private islands were detected for the Hungarian sample on ECA8, ECA14, ECA23 and ECA30 and for the Austrian sample on ECA3, ECA5, ECA8, ECA16 and ECA20. Within the Austrian Lipizzans an overlapping 953.3kb long homozygous region on ECA16:37.87-38.82 Mb was found, that contained among 29 annotated genes the two genes ECATH-3 and ECATH-2. Both were described as functional genes encoding Cathelicidin-derived antimicrobial peptides, involved in the peptidebased host defence of neutrophils and epithelia in horses [21,22]. GO analysis for the Austrian Lipizzan gene list highlighted the term GO:0006950~response to stress for the gene UCN2 (urocortin 2, a member of the corticotropin-releasing hormone family) on ECA16 and the gene HSP90AB1 (Equus caballus heat shock protein 90 alpha family class B member 1) on ECA20, whereas the latter was assigned to four Equus caballus species by Vidale et al. [23]. HSP90AB1 is considered to play a role in the age-dependent metabolism of chondrocytes in horses [24].
All subpopulations shared a homozygous region on ECA11, nearly identical in position and length (mean length 676.4 kb) at 24.13-24.81 Mb containing the annotated genes COPZ2, CBX1, SNX11, HOXB1, HOXB2, HOXB3, HOXB5, HOXB6, HOXB7, HOXB8, HOXB13, TTLL6 (and CALCOCO2, additional in the Hungarian sample). For the Croatian/Hungarian Lipizzans an additional island on ECA4:58.10-58.60 Mb was confirmed containing the homeobox-A cluster. Enrichment analysis highlighted for all samples highly significant levels (p < 0.001) for the homeobox-B cluster, mostly related to biological processes involved in embryonic morphogenesis and anterior/posterior specification (GO:0009952, GO: 0048704, GO:0009953, GO:0009954). Significant levels for GO terms based upon the HOXA cluster (GO:006 0065~uterus development, GO:0009954~proximal/distal pattern formation, GO:0007338~single fertilization, GO:00 08584~male gonad development) were exclusively documented for Croatian and Hungarian Lipizzans. The island on ECA11 including the HOXB cluster and highly significant GO terms related to morphogenesis were recently also documented in the Posavina horse by Grilz-Seger et al. [16]. Zhang et al. [25] listed HOXB13 in a recent selection signature study in horses. HOX genes play a fundamental role for morphological diversity in animals and for the control of axial morphology along the anterior-posterior body axis. As these gene clusters are directly assigned to axial body segments in numerical order, they can control single segments during embryogenesis concerning their position, segmentation and further differentiation [26].
On ECA11:31.01-31.94 another common homozygous region was found in all samples, whereas the frequency in Hungarian Lipizzans was with 47% slightly below the threshold. This ROH island, nearly identical in length (926.9/910.2 kb) in the Austrian and Croatian Lipizzans, harboured the genes C11H117orf67, DGKE, COIL, SCPEP1, AKAP1, MSI2 and was partially overlapping with two smaller islands of the Slovakian population. The Hungarian sample had a 244.4 long overlapping homozygous region within the Musashi RNA binding protein 2 (MSI2), which is involved in stem cell proliferation [27], haematopoiesis and can promote aggressive myeloid leukaemia in humans [28]. Four of these genes (COIL, SCPEP1, AKAP1, MSI2) were listed among other 121 genes that underwent positive selection during the domestication process of the horse [29].
In general the majority of ROH islands were medium sized (between 100 and 900 kb) to small sized (8-100 kb) according to the classification of Pemperton et al. [30]. Only in the Croatian sample one island on ECA14 was longer than 1 Mb. Most of the short ROH islands did not exceed the length of the contained gene, e.g. the island on ECA3:118.66-118.76 Mb, containing the two genes UV-stimulated scaffold protein A (UVSSA) and macrophage erythroblast attacher (MAEA). UVSSA is involved into the repair process of DNA damaged by ultraviolet (UV) sun rays. These two genes were directly followed by a 83 kb long island containing the gene C-terminal-binding protein 1 (CTBP1), a transcriptional co-repressor. Pemperton et al. [30] emphasized the investigation of short ROHs (< 100 kb), as they reflect homozygosity for ancient haplotypes contributing to local linkage disequilibrium (LD) patterns. Intermediary ROHs longer than 100 kb up to 2.5 Mb are deemed to be the result from background relatedness in populations with limited census. Furthermore, Pemperton et al. [30] stated that by the commonly applied windows approach as implemented in PLINK [31], ancient haplotypes and small ROHs due to high recombination rate might remain undetected. With the 500 kb applied boundary we were able to detect very small ROH islands, but SNP extraction revealed, that the frequencies of animals sharing those islands were underestimated, exemplarily demonstrated for the genes MAEA, UVSSA, CTBP1 and the STX17 locus.
Grey coat color, caused by a 4.6 kb duplication in intron 6 of Syntaxin 17 (STX17), is an essential part of the breeding objective in Lipizzan horses and it has been under selection since the past 150 years. The STX17 mutation is embedded into 352 kb long haplotype [13,32].The frequency of the identified ROH island (ECA25:6,39-6,79) harbouring the genes STX17, NR4A3, ERP44 and INVS, was in concordance with genotype frequencies cited in Pielberg et al. [13], and ranged from 51.1% (Croatia) to 71.2% (Austria) among the respective subpopulations. In the Austrian breeding program only grey horses are used for reproduction, whilst the Slovakian, Hungarian and Croatian stud farms also permit solid colored horses. As a consequence the proportion of homozygous grey G/G horses is higher in the Austrian population.
The allele status of STX17 influences the progression of greying by age, grade of speckling (pigmentation among grey hair), grade of vitiligo and incidences and grade of melanoma [10][11][12][13]. The incidence of melanoma in Lipizzans was reported to affect 80% of the animals older than 15 years [10,13]. Mostly melanomas in Lipizzans are benign and do not metastasize [33]. A previous study revealed, that the genes STX17, NR4A3, ERP44 and INVS, which are located in the aforementioned ROH island also were expressed in melanoma tissues [13]. Pielberg et al. [13] observed an over-expression of NR4A3 (nuclear receptor subfamily 4, group A, member 3), which is influenced by the dominant G-allele. The underlying molecular processes are not clearly understood [34]. Seltenhammer et al. [33] investigated factors retarding metastasis in affected grey Lipizzan horses and revealed common features with blue nevi in humans and identified potential markers for equine melanoma (S-100, PCNA, HMB-45, Ki-67, T-311, CD44). We identified within the entire Lipizzan sample an island on ECA14 containing the genes SPRY4 (Sprouty4) and NDFI1 (Nedd4-family interacting protein 1). SPYR4 is a member of the SPYR proteins, which are negative regulators of growth factor signaling pathways. Thus, SPRY4 inhibits transformed cell growth, migration and invasion, and epithelial-mesenchymal transition [35]. Shaverdashvili et al. [36] provided evidence in human melanoma for an inverse relationship between the expression of SPYR4 and Membrane-type 1 Matrix Metalloproteinase (MT1-MMP), one important driver of melanoma metastasis. Inhibition of MT1-MMP resulted in increased expression of SPYR4 and was correlated with the survival of melanoma patients in humans. NDFIP1 also acts as a tumour suppressor by repressing cell proliferation and was reported to be down-regulated in human uveal melamoma [37]. In a previous study Jiang et al. [34] investigated the ERK pathway and studied the effects of BRAF, RAS, GNAQ, GNA11, KIT on melanoma development in grey horses, whereas an association between STX17 and an activation of the ERK pathway was shown, but no activation of the aforementioned genes were observed. In the Austrian Lipizzan sample we identified an island containing the gene HSP90AB1. Metri et al. [38] identified HSP90AB1 as a main discriminator between metastatic and primary melanoma in humans, which also was correlated with survival rate of melanoma patients. In our data GO analysis revealed an enrichment of HSP90AB1 in the terms GO:0069550-reponse to stress and GO:0071353-cellular response to interleukin-4.

Conclusions
The enhanced high-resolution Netview approach results in a generally understandable visualisation, which is able to simultaneously present fundamental metrics of genetic diversity, coancestry and relatedness on population and individual level. This approach thus provides a high practical impact for conventional and conservation breeding management. The analyses of ROH islands based upon the conventional ROH size boundary of > = 500 kb have proven to be a valuable approach to screen for signatures of selection, also able to detect short homozygous overlapping regions. In order to avoid an underestimation of the frequencies of very short ROHs the underlying size boundary should be adapted in regions of special interest. The investigation of subpopulations within one single breed allows a fine-scale analysis and we demonstrated that slightly different breeding objectives have a high impact on genomic regions and gene content. ROH island analyses identified several annotated genes (SPRY4, NDFIP1, HSP90AB1 and IMPDH2), which may play a role for further studies on equine melanoma.

Ethics statement
This study was discussed and approved by the institutional Commission for Ethics and Animal Welfare, University of Veterinary Medicine, Vienna, protocol number: ETK-06/ 05/2015, in accordance with GSP guidelines and national legislation. The Lipizzan state stud farms Piber (Austria), Topol'čianky (Slovakia), Lipik, Đakovo (Croatia) and Szilvasvárad (Hungary) granted the permission to take hair samples from their horses.

Sampling
The horses included in this study originated from five European state stud farms, which are located in four different countries. The Austrian and the Slovakian samples represent the entire breeding populations of the stud farm of Piber (Austria, 254 horses) and Topol'čianky (Slovakia, 55 horses). The Croatian Lipizzans were sampled from the state stud farms Lipik (8 horses) and Đakovo (37 horses), whilst 23 Lipizzan horses were collected from the Hungarian state stud farm Szilvasvárad.
The hair samples were taken between the years 2013 and 2017 and the selected horses were born between 1989 and 2014. For Austria and Slovakia we were able to collect the whole stud farm populations. The breeding stock of the Croatian stud farms comprises together approximately 120 breeding animals, the number in the Hungarian state stud farm is on a similar level.

SNP genotyping
The single nucleotide polymorphism (SNP) genotypes for the 377 horses were determined using the Affymetrix Axiom™ Equine genotyping array [14]. The chromosomal position of the SNPs was derived from EquCab2.0 reference genome. We excluded SNPs positioned on the sex chromosomes (X: 28,017 SNPs and Y: 1 SNP), SNPs without known chromosomal position (30,864 SNPs) and SNPs with more than 10% missing genotypes, resulting in a total of 611,914 SNPs for each horse. For the analysis of population structure and genetic diversity we excluded SNPs with a minor allele frequency (MAF) less than 0.01, resulting in 492,719 SNPs per horse.

Genetic diversity and ROH analysis
In order to illustrate the population structure we applied Principal Component Analysis (PCA) on basis of the genetic relationship matrix (G) with pairwise identities by state (IBS) between horses as provided by PLINK v.1.7 [31]. The PCA was performed using R platform, whilst pairwise F ST values were calculated with the R package Geneland (https://rdrr.io/cran/Geneland).
Population clustering and admixture analysis were determined using the program Admixture 1.23 [39]. We ran Admixture for 100 iterations increasing K from 2 to 10. 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 variation between runs. Cross validation (CV) error estimation for each K was performed to determine the optimal number of clusters. Admixture results were visualized with the program Distruct 1.1 [40] and integrated in the high-resolution network following the procedure of Druml et al. [15].
ROH were determined with an overlapping window approach implemented in PLINK v1.7 [31] based on the settings: minimum SNP density = one SNP per 50 kb; maximum gap length = 100 kb; minimum length of homozygous segment > 500 kb (including more than 80 homozygous SNPs); one heterozygote and two missing genotypes were permitted within each segment.
The total number of ROH (N ROH ), average length of ROH (L ROH ) and the sum of all ROH segments (S ROH ) were summarized according to respective subpopulations. To analyze ROH length classes, ROHs were divided into following seven length categories: 0.5-1, 1-2, > 2-4, > 4-6, > 6-8, > 8-10 and > 10 Mb. The genomic inbreeding coefficients (F ROH ) were calculated by use of following formular: ,whereas the length of the autosomal genome (L AUTO ), was set to 2243 GB. Statistical analyses and graphical representations of specific genomic regions and distribution of ROH segments were performed using the R-package detectROHs (www.r-project.org)and SAS v.9.1. [41].

High-resolution population networks
In order to ascertain the high-resolution population structure of Lipizzan horses, we performed a network visualization based upon the aforementioned IBS derived relationship matrix (G), following the description of Druml et al. [15]. The different components involved in the so-called NetView approach are described in detail by Neuditschko et al. [42] and Steining et al. [43]. Briefly, we computed genetic distances by subtracting pairwise relationships from 1 and applied the algorithm in its default setting (number of k nearest neighbors k-NN = 10). To illustrate the genetic relatedness between neighboring horses, we associated the thickness of edges (connecting lines) with the proportion of the genetic distance, whilst thicker edges corresponding to lower genetic distances. To identify highly inbred and outcrossed horses within the respective population networks we scaled the node size of each horse based on the individual S ROH . The node colors of each horse represent the individual level of admixture at the selected number of K-clusters, whilst the node border colors illustrate the origin of the horses.

Gene characterization in ROH islands
The distribution of ROH segments across the genome was visualized using the R-package detect ROHs (www.r-project.org). Putative ROH islands were determined based upon overlapping homozygous regions, which were shared by more than 50% of the horses. ROH islands were calculated for the entire data sample and respective subpopulations. Genes located in ROH regions were identified with the map viewer of the equine ensemble database EquCab2.0 (www.ensembl.org) and Gene Ontology (GO) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways of annotated genes were analyzed using the open source Database for Annotation, Visualization, and Integrated Discovery (DAVID) v6.8 package (https://david.ncifcrf.gov) using the Equus caballus annotation file as background [44].