Genomic inbreeding trends, influential sire lines and selection in the global Thoroughbred horse population

The Thoroughbred horse is a highly valued domestic animal population under strong selection for athletic phenotypes. Here we present a high resolution genomics-based analysis of inbreeding in the population that may form the basis for evidence-based discussion amid concerns in the breeding industry over the increasing use of small numbers of popular sire lines, which may accelerate a loss of genetic diversity. In the most comprehensive globally representative sample of Thoroughbreds to-date (n = 10,118), including prominent stallions (n = 305) from the major bloodstock regions of the world, we show using pan-genomic SNP genotypes that there has been a highly significant decline in global genetic diversity during the last five decades (FIS R2 = 0.942, P = 2.19 × 10−13; FROH R2 = 0.88, P = 1.81 × 10−10) that has likely been influenced by the use of popular sire lines. Estimates of effective population size in the global and regional populations indicate that there is some level of regional variation that may be exploited to improve global genetic diversity. Inbreeding is often a consequence of selection, which in managed animal populations tends to be driven by preferences for cultural, aesthetic or economically advantageous phenotypes. Using a composite selection signals approach, we show that centuries of selection for favourable athletic traits among Thoroughbreds acts on genes with functions in behaviour, musculoskeletal conformation and metabolism. As well as classical selective sweeps at core loci, polygenic adaptation for functional modalities in cardiovascular signalling, organismal growth and development, cellular stress and injury, metabolic pathways and neurotransmitters and other nervous system signalling has shaped the Thoroughbred athletic phenotype. Our results demonstrate that genomics-based approaches to identify genetic outcrosses will add valuable objectivity to augment traditional methods of stallion selection and that genomics-based methods will be beneficial to actively monitor the population to address the marked inbreeding trend.

the most comprehensive genetic analysis performed in this population (n = 10,118) we show a striking temporal increase in inbreeding and regional variance across the global Thoroughbred population during the last five decades. Individual inbreeding coefficients (F IS ) estimated using 9,212 pruned SNPs and runs of homozygosity (F ROH ) (minimum length 1 Mb) characterised using an unpruned set of 46,478 SNPs revealed a highly significant increase in inbreeding over time (F IS : R 2 = 0.942, P = 2.19 × 10 −13 ; F ROH : R 2 = 0.88, P = 1.81 × 10 −10 ) (Fig. 2) with the greatest rate of change observed since the 2000s (S12 Figure). A linear regression model was used to test for significance and directionality of change in annual mean inbreeding with respect to year of birth (S12 Figure). Inbreeding estimates were determined using two methods: individual inbreeding coefficients (F IS ) measure observed versus expected genetic diversity in an individual in a population. F ROH estimates the proportion of the genome covered by runs of homozygosity. The values are correlated but the unit of measurement is different so the absolute values cannot be directly compared. Similar trends for both measures were observed within each geographic region (S13-S16 Figures). Results from Student's t-tests indicated that inbreeding estimates were higher in EUR compared to ANZ (F IS , P = 1.9 × 10 −21 ; F ROH , P = 0.043) and NAM (F IS , P < 3.1 × 10 −19 ; F ROH , P = 0.00059). A temporal reduction of Thoroughbred genetic diversity has previously been reported among much smaller sample sets from single geographic regions 8,9 . Expansion of the timeline to include the most recent decade using a sample size more than 20-fold larger indicates that despite industry cognizance of inbreeding and previous cautions 8 , there has been no arrest in the rate of increase in inbreeding and it is a global, population-wide phenomenon.
Breeding practices that promote inbreeding have not resulted in a population of faster horses 4,10 and our results, generated for the first time using a large cohort of globally representative genotypes, corroborate this. When evaluating Timeform handicap ratings for EUR horses, no association was found between F IS and racing performance for horses born between 1996-2013 (n = 1,886, r 2 = 0.034, P = 0.780) and 2009-2013 (n = 1,065, r 2 = −0.014, P = 0.65). Purposeful inbreeding attempts to duplicate favourable gene variants that are selected at each successive generation, but homozygosity arising from inbreeding tends to be associated with decreased trait values since the proportionate trait gain arising from beneficial mutations is generally limited by pleiotropy 11 . Also, among domesticates, there is an increased mutational load and a higher proportion of deleterious alleles in regions under selection [12][13][14][15] . Therefore, there may be a genetic cost of inbreeding, which in dogs has negatively influenced intra-breed genetic diversity and increased the frequency of disease-causing alleles 15,16 and in cattle, is associated with decreases in milk, fat and protein yields and negatively impacts reproductive traits 17 .
While inbreeding depression generally negatively impacts health and fertility through the accumulation of deleterious alleles, the 'tipping point' at which there is an irreversible accumulation of unfavourable mutations is not currently known. In horses mutational load and purging of deleterious alleles has been assessed in whole genome sequence data and revealed a significant relationship between inbreeding and mutational load 18 . Year of Birth Surprisingly, according to a recent study 18 , Thoroughbreds appear to have a lower than expected mutational load, suggested to be due to effective purging through negative selection on phenotypes. This may be facilitated in the Thoroughbred in practice by the unusually large census population size relative to the effective population size, with a high proportion of horses that do not ever race 19 . However, the results from that study are likely not representative of the genetics of the current breeding population in the major bloodstock regions of the world; the sample cohort of Thoroughbreds examined was small (n = 19) and all the horses were registered Thoroughbreds in Korea. For example, ROH identified in that study were up to 11 Mb, whereas we have detected here ROH up to 40 Mb (S1 Table); increasing ROH increases the likelihood of deleterious alleles being exposed. It is interesting to consider also that many of the performance limiting genetic diseases in the Thoroughbred do not generally negatively impact on suitability for breeding; some diseases, with known heritable components, are successfully managed by surgery (osteochondrosis desicans, recurrent laryngeal neuropathy), nutritional and exercise management (recurrent exertional rhabdomyolysis) and medication (exercise induced pulmonary haemorrhage). This facilitates retention of risk alleles in the population and enhances the potential for rapid proliferation of risk alleles if they are carried by successful stallions. Furthermore, since single gene variants have not yet been identified for these diseases, it is likely that they result from polygenic additive genetic variation that may not be easily exposed and purged by negative selection on individual ROH. In order to fully understand the consequences of inbreeding in the Thoroughbred, inbreeding measures should be regressed on large population-scale phenotypes and the extent of mutational load should be determined in a large cohort of representative samples.
There are currently concerns about the impact of the increasing use of small numbers of stallions in the breeding population on inbreeding and population viability. While the census population size of the Thoroughbred is large, a more appropriate method to assess population health is to estimate the effective population size (N e ), an estimate of the size of an idealised population that is representative of the genetic diversity in the actual population [20][21][22] . To better understand the impact of inbreeding on diversity in the population we estimated N e for the global population and regional populations using the subset of horses born 2013-2017 (n = 3,341) to represent the current breeding population. The global population had N e = 330 and among regional cohorts N e was highest in NAM (N e = 226) and lowest in SAF (N e = 93); also ANZ (N e = 197), EUR (N e = 198). The observation that N e was higher when all regions were considered together indicates that there is opportunity to exploit regional variation to improve genetic diversity in the global population. This may be achieved in practice, for example, by selecting stallions that are genetically distant from mares, or more broadly by the permanent movement of stallions to regions that have a genetically diverse population of mares. A comparison of N e with census population sizes (N c , included here as N c = 100,000 for NAM, EUR, ANZ and SAF and N c = 500,000 global) indicates that the genetic variation in the population deviates substantially from expectations. N e /N c was <0.2% in all regional cohorts and was <0.1% when the population was considered as a whole. In most wild species N e is 10-50% of the census population size and among domesticates N e /N c has a median of 3% 23 . A limitation to interpretation of the results from this study is the potential bias in the samples used, many of which came from the same breeding farms. Also, in some cases, the sample sizes are relatively small (e.g. SAF) and may not be representative of the entire population. Nonetheless, these results highlight the need for a systematic evaluation of genetic diversity that may be applied for longitudinal monitoring of the Thoroughbred population.
Selection for the thoroughbred phenotype. As well as the effect of close relationships among breeding individuals, inbreeding is a consequence of selection for favourable traits. Here, to systematically identify genes that have been targets of selection for the Thoroughbred athletic phenotype, and may be most impacted by inbreeding, we used the composite selection signal (CSS) approach, a weighted measure that combines the signals from identification of highly differentiated loci (F ST ), increase in selected allele frequency (ΔSAF) and cross-population extended haplotype homozygosity (XP-EHH) tests 24,25 , to analyse 48,896 pan-genomic SNPs genotyped in elite Thoroughbreds (n = 110) and non-Thoroughbreds (n = 84) representing putative founder populations. We identified 15 significant candidate selected genomic regions (S2 Table, S17 Figure), defined as clusters of ≥5 SNPs among the top 1% of the smoothed CSS statistic result (−log 10 P), seven of which overlapped with previously reported genomic regions with evidence for selection in Thoroughbreds 26,27 (Table 1). As a high prevalence of runs of homozygosity (ROH) can inform on selection, the top 1,000 SNPs ranked by the percentage of individuals with a certain SNP located within a ROH were extracted for ROH >1 Mb and >5 Mb (S3-S4 Tables, S18 Figure). Among the CSS regions, there was substantial overlap with regions with a high prevalence of ROH.
We interrogated the CSS peaks as well as flanking regions (±0.5 Mb) for candidate genes that may contribute to the Thoroughbred phenotype and identified genes with functions in behaviour, musculoskeletal, cardiac and respiratory function, conformation and metabolism ( Table 1). Given that positive selection at a specific genomic locus tends to reduce ('sweep') variation across a larger region, it can be difficult to identify the gene under selection. Here, we identified plausible candidate genes driving selection based on the location of the highest-ranking SNPs in the selected region and by reviewing known biological functions of genes that we hypothesise may be selected for the Thoroughbred phenotype. Seven of the top 10 ranked SNPs were located in the top ranked region on ECA1 and the top three SNPs spanned 95 kb closest to the ZW10 interacting kinetochore protein gene (ZWINT). ZWINT (also known as SIP30) is abundantly expressed in the brain, modulates neurotransmitter release and functions in the mediation of peripheral nerve injury-induced neuropathic pain [28][29][30] . In rodents ZWINT influences pain-evoked emotional responses 31 and since human athletes have been reported to have a greater ability to tolerate pain than normally active controls 32 , it is intriguing to speculate ZWINT may be involved in the equine response to exercise-induced pain. The region containing ZWINT was also identified as having reduced genetic diversity in Japanese Thoroughbreds compared to other breeds 27 . There were particularly long ROH in this region in the current study; 31 of the 305 stallions had ROH >16 Mb between ECA1: 37.8-76.3 Mb (S1 Table); the longest ROH spanned almost 40 Mb. (2020) 10:466 | https://doi.org/10.1038/s41598-019-57389-5 www.nature.com/scientificreports www.nature.com/scientificreports/ We identified other neurological/behaviour associated genes in the selected regions including the follistatin like protein 4 gene (FSTL4), the abnormal spindle microtubule assembly gene (ASPM) and the X-prolyl aminopeptidase 1 gene (XPNPEP1). Knockdown of FSTL4 in mice results in the extinction of inhibitory avoidance memory indicating its involvement in synaptic plasticity and memory formation. Interestingly it functions by directly interacting with the exercise-induced brain derived neurotrophic factor (BDNF) 33,34 which decreases in response to chronic stress 35 . ASPM is a major determinant of the size of the cerebral cortex in mammals 36,37 that plays a key role in memory, attention, perception and awareness. Knockdown of XPNPEP1, which encodes an important downstream regulator of the stress response 38 , in mice results in enhanced locomotor activity and impaired contextual fear memory 39 . An emerging theme in equine transcriptomics and genomics research suggests a link between the exercise phenotype and behavioural plasticity. For example, in the skeletal muscle transcriptome response to exercise training, neurological processes were the most significantly over-represented gene ontology (GO) terms, with the top three ranked GO terms being Neurological system process (P = 4.85 × 10 −27 ), Cognition (P = 1.92 × 10 −22 ) and Sensory perception (P = 4.21 × 10 −21 ) 40 . Furthermore, in genome-wide association (GWA) studies genes involved in behavioural plasticity are the most strongly associated with economically important traits in racing Thoroughbreds: precocity (early adaptation to racing) 41 and the likelihood of racing 19 . For Thoroughbred horses, behavioural plasticity enables adaptation to the rigours of an intense exercise training programme in an unnatural environment, with considerable variation in the abilities of horses raised in the same environment to adapt to stress. The presence of these genes in genomic regions under selection in the Thoroughbred population supports human-mediated adaptation of the Thoroughbred towards heightened awareness and the ability to learn and adapt to stress.
Genes for aesthetic physical phenotypes including height, stature, coat and plumage colour are commonly encountered in genomic regions under selection in domestic animal populations since they are easily identifiable [42][43][44][45] . Conformation characteristics are among the most discernible traits in horses and are the principal observable traits on which Thoroughbreds are selected. Here, the selection signal on ECA21 centred on the iroquois homeobox 1 gene (IRX1) and the ADAM metallopeptidase with thrombospondin type 1 motif 16 gene (ADAMTS16), a locus that has been shown to be associated with hip geometry in humans 46 . While ADAMTS16 has been identified among selection signals for high altitude adaptation in pigs 47 , IRX1 is associated with bone mineral density in humans 46 and influences chondrocyte differentiation and may therefore contribute to joint www.nature.com/scientificreports www.nature.com/scientificreports/ flexibility 48 . In horses, genes with functions in bone mineral density have been associated with measurements of joint angles 49 . Since the angle of the pelvis is a major determinant of physical conformation in Thoroughbreds and is associated with injury and performance 50 , we hypothesise that the selection signature on ECA21 reflects the evolution of the Thoroughbred conformation phenotype.
Other genes relating to musculoskeletal form and function were identified among selected regions. For example, the deleted in lymphocytic leukemia 7 gene (DLEU7) is associated with growth traits in chickens 51 and overgrowth in humans 52 and may therefore contribute to stature in Thoroughbreds. In bone physiology, the von Willebrand factor C domain containing 2 gene (VWC2), a bone morphogenic protein, promotes bone formation, regeneration and healing 53,54 . In the context of Thoroughbred musculature, the multiple EGF like domains 10 gene (MEGF10) controls muscle cell proliferation and is involved in the regulation of myogenesis 55,56 and the collagen type VI alpha 3 chain gene (COL6A3) plays a major role in the maintenance of strength of muscle and connective tissue 57 . COL6A3 is one of three genes encoding components of collagen VI, which in the horse is expressed in developing cartilage 58 . Also, collagen VI disruption in horses has been associated with osteochondrosis, a common developmental orthopaedic disease with a major economic impact in the Thoroughbred industry 59 . COL6A3 may also be relevant to the muscle metabolism phenotype since its expression in adipocytes is associated with insulin resistance and obesity 60,61 . Mutations in the diacylglycerol O-acyltransferase 2 gene (DGAT2) mutations cause Charcot-Marie Tooth disease in humans 62 and the DGAT2 protein also functions in lipid metabolism 63-67 . Polygenic adaptation in the Thoroughbred. There are likely numerous endophenotypes on which selection acts to generate the athletic phenotype. Adaptation driven by selective sweeps at a number of key genomic loci is likely to be important; however, modest allele frequency changes at multiple loci are also expected to occur as a consequence of polygenic adaptation 68,69 . Therefore, in addition to 'core' genes that are critical to the phenotypic outcome, highly granular additive genetic variation-essentially encompassing the entire genomecombined with epistatic interactions differing across cell types, may largely shape the phenotype 70 . Considering this, we performed an enrichment analysis to identify functional processes over-represented among the set of 387 of the 462 genes contained within the putative selection regions that mapped to the Ingenuity Pathway Analysis (IPA) database. The top canonical pathway was Airway inflammation in asthma (S5 Table). Manual curation of the top 50 canonical pathways points to a process of polygenic adaptation among functional modules in cardiovascular signalling, cellular growth, proliferation and development, cellular immune response, cellular stress and injury, metabolic pathways (fatty acid and lipid degradation/biosynthesis), neurotransmitters and other nervous system signalling and organismal growth and development (S6 Table). Our results support the development of the Thoroughbred phenotype via a contribution from major allele frequency shifts at 'core' genes contributing to behavioural, metabolic and conformation traits and genome-wide changes in functional modules that shape a range of exercise-relevant physiological adaptations.
Concluding remarks. We report here a highly significant increase in inbreeding in the global Thoroughbred population during the last five decades, which is unlikely to be halted due to current breeding practices. Inbreeding results in mutational load in populations that may negatively impact on population viability. 'Genetic rescue' of highly inbred populations may be possible by the introduction of genetically diverse individuals 71 ; however, rescuing genetic diversity in the Thoroughbred will be challenging due to the limitations of a closed stud book. Furthermore, the population has a small effective population size (N e ) and a limited numbers of stallions have had a disproportionate influence on the genetic composition of the Thoroughbred; 97% of pedigrees of the horses included here feature the ancestral sire, Northern Dancer (1961) and 35% and 55% of pedigrees in EUR and ANZ contain Sadler's Wells (1981) and Danehill (1986), respectively (S1 Text).
Pedigree data can be useful to illustrate broad trends in breeding practices, but since pedigree-based estimates of inbreeding and relatedness have poor correlations with estimates using genomic methods 27,72-75 relying on pedigree alone for outcrossing is likely to be inefficient in reversing the trends observed here (S1 Text, S19 Figure). Directives to prevent over-production from popular sire lines and the global movement of stallions that are distinct from the local population of mares may act to maintain and increase genetic diversity in the population. However, given the limited diversity in current Thoroughbred pedigrees, genomics-based measures using high-density genome-wide SNP information and a large reference population are likely to offer the best opportunity to slow and reverse the potential effects of inbreeding. The introduction of an industry-mediated longitudinal programme of genomics-based monitoring of inbreeding and the implementation of guidelines and strategies for genome-enabled breeding that are comparable to methods used in other domestic species, will contribute to promoting economic gain and safeguarding the future of the breed.

Methods
Ethics statement. Samples from animals used in this study were collected by owners and submitted to Plusvtial Ltd. for commercial genetic testing. Consent for use of samples in research was obtained during the sample submission process and methods were carried out in accordance with the agreement. No experimental procedure was performed on live animals.
Sampling and population assignment. The following data was collected for n = 10,118 Thoroughbred (TB) horses: horse name, sire, dam, sex, year of birth and country of birth (S7 Table). Based on country of birth, horses were assigned to the following geographic regions: Europe/Middle East (EUR), Australasia (ANZ), North America (NAM), South Africa (SA) and Japan (JAP). An additional set of n = 84 horses from four breeds that were chosen to represent putative TB founder populations included n = 20 Akhal Teke (AKTK), n = 24 Arabian (ARR), n = 23 Moroccan Barb (MOR) and n = 17 Connemara pony (CMP). The ARR and AKTK sample data used were obtained from publicly available equine Illumina SNP50 Beadchip genotype data 3 .
Assembly of comparative SNP data set, quality control and filtering of SNPs. DNA was isolated from blood or hair samples and genotyped using the Illumina EquineSNP50 BeadChip (SNP50), the Illumina EquineSNP70 BeadChip (SNP70) or the Affymetrix Axiom ™ Equine 670 K SNP genotyping array (SNP670). Only animals and SNPs with a genotyping rate >95% were included with a minor allele frequency (MAF) threshold >0.05 applied. A set of 48,896 autosomal SNPs derived originally from the SNP50 and SNP70 arrays was used for the analysis. This SNP set was extracted from the genotype data from each of the three platforms. 1,821 SNPs were not present on the SNP670 array. SNPs that failed quality control in <5% of samples or were not present on one of the array platforms were imputed using the software program BEAGLE (version 3.3.2) 76,77 . For 10 horses genotyped using both the SNP50 and SNP70 arrays and 10 different horses separately genotyped using the SNP70 and SNP670 array post-imputation concordance was greater than 99%.The TB dataset (n = 10,118) comprising of the set of 48,896 SNPs was pruned using the PLINK software suite V1.9 (http://pngu.mgh.harvard.edu/purcell/plink/) 78 -indep function with the following parameters: a five-step sliding window size of 50 with a VIF threshold of 50 where the VIF is 1/(1 − R 2 ). This pruned set of 9,212 SNPs was used for principal component analysis (PCA) within the Thoroughbred population and for the calculation of individual inbreeding coefficients (F) as outlined below.
The Thoroughbred population was compared to founder populations using PCA and composite selection signature (CSS) analysis. We randomly selected n = 229 elite horses (CPI > 5, i.e. earned more than five times the average; CPI is a class performance index defined by the American Jockey Club) from the TB dataset (n = 10,118). Genomic relationships among all horses, Thoroughbreds and non-Thoroughbreds, were estimated using autosomal identity by descent (pi-hat) values in PLINK v1.9 79 . After removing individual horses with pi-hat > 0.25, 150 horses (n = 23 MOR, n = 17 CMP, and n = 110 TBE) remained. Then, the newly genotyped data was merged with publicly available data for AKTK (n = 20) and ARR (n = 24). Each population was pre-processed separately to include only individuals and SNPs with a genotyping rate > 95% and with MAF > 0.05. Following breed specific quality control a final round of quality control (MAF > 0.01 and genotyping rate >0.95) was applied to the combined data set of 194 horses. This resulted in 31,722 SNPs remaining for the mixed breed PCA and CSS analyses. The datasets created are summarised in Table 2. The sex, year of birth and region of origin of the Thoroughbred samples are provided in S7 Table. Stallion genotype reconstruction. To increase the representation of prominent stallions genotypes were reconstructed for 127 horses. Of these 43 had previously been genotyped and were used to assess the accuracy of the genotype reconstruction (S1 Text). The genotypes (~46,000 autosomal SNPs) were reconstructed for horses where genotypes were available for n ≥ 20 progeny of an individual stallion (S8 Table). An adaptation of the method described by Gomez-Raya 80 was used to infer genotypes. The main difference here was that population genotype frequencies were included allowing accurate reconstruction with just 20 offspring. If all offspring do not share one allele at a locus the sire must be heterozygous at the locus. If all offspring do share one allele at a locus the sire may be either heterozygous at the locus or homozygous for the shared allele. However, the probability of the sire being homozygous or heterozygous can be calculated based on the proportion of the different genotypes in the offspring and the observed proportion of each genotype in the population. The chi-square statistic was calculated using observed and expected offspring allele frequencies for each of the three possible sire genotypes. In order to compare the three sire genotype possibilities and assign a relative probability to each, the chi-square test statistics were first converted to likelihood scores. To do this the density of a chi-square distribution was taken at each of the three points (the density at each point represents the likelihood of the sire having that genotype, given the observed offspring genotypes). The density values were then divided by the sum of the three densities to normalize. Normalizing these three likelihoods by their sum gave the relative likelihood of each sire genotype and the genotype with the maximum likelihood was assigned as the sire genotype.
The option outliersigmathresh is used to identify samples which exceed a defined number of standard deviations along one of the top principal components and classify as outliers. The threshold for this was increased  www.nature.com/scientificreports www.nature.com/scientificreports/ from the default threshold of six to a threshold of 10 to ensure that samples from the Founder populations were not mis-classified as outliers. All other parameters were set to default values. Principal components (PCs) were plotted with data points colour-coded based on each horse's geographic region of origin. A series of analyses were performed to investigate Thoroughbred population sub-structure using the Thoroughbred and Founders, the Pruned Thoroughbred and the Pruned Sire set as described above and summarised in Table 2.
Inbreeding estimates. Using the Pruned Thoroughbred Set individual inbreeding coefficients (F IS ) for each horse were calculated in PLINK based upon reduction in heterozygosity relative to Hardy-Weinberg expectation (thehet option in PLINK). Genomic inbreeding was also evaluated in the Unpruned Thoroughbred Set by identifying runs of homozygosity (ROH) with the-homozyg command. ROH were defined as tracts of homozygous genotypes that were >1 Mb in length identified for one SNP per 1000 kb on average and two consecutive SNPs <1000 kb apart. No more than two missing genotypes and one heterozygous genotype were allowed. The following parameters were set:-homozyg;-homozyg-kb 1000;-homozyg-snp 30;-homozyg-gap 1000;homozyg-window-het 1;-homozyg-window-snp 30;-homozyg-density 1000;-homozyg-window-missing 2 to identify runs of homozygosity spanning at least 1 Mb. The threshold for the number of SNP per RoH was set at 30 in order to identify shorter RoH of 1-2 Mb. As over 80% of 50 SNP windows span a region of greater than 2 mb these shorter RoH could not be detected using a minimum number of 50 SNPs per RoH. The analysis was also run with the parameter homozyg-kb increased to 5000 to identify runs of homozygosity of 5 Mb or greater reflective of more recent inbreeding 82 (Keller 2011). First, the individual sum of total ROH per animal was calculated. The F ROH statistic proposed by McQuillan et al. 83 was then calculated, whereby the total length of ROH covering an individual animal's genome (L ROH ) is divided by the length of the autosomal genome (L AUTO ); F ROH = L ROH / L AUTO . Here, consistent with other equine studies 84,85 we used the length of the equine autosomal genome (assembly EquCab 2) 86 as 2,242,960 kb (www.ncbi.nlm.nih.gov/genome/145?genome_assembly_id=22878).
The annual mean, SE, SD and CI were calculated based on the horse's year of birth for both measures of inbreeding (S16-S17 Table). Following a preliminary analysis, the pre-1996 samples were excluded from further analysis as sample distribution pre-1996 was low. The post-1995 data is of the most interest as this coincides with the introduction of "big books" for stallions; i.e. large numbers of mares bred. The shuttling of stallions to Australia from Europe for dual hemisphere breeding seasons peaked in 2001. A summary of the year of birth, region of birth and sex of the horses used in these analyses is provided in S7 Table. A pair-wise Student's t-Test was used to compare measures of inbreeding across the main geographic regions represented in the data set; Europe/ Middle East (EUR), Australasia (ANZ) and North America (NAM).
A linear model was used to assess the relationship between year of birth and inbreeding for the global population and for each of the geographic regions. Mean annual inbreeding values were used in the regression model (S9 Table, S10  Table). Pearson correlation was run to assess the relationship between inbreeding and racing performance defined by handicap (Timeform) rating. Significance was calculated by testing the null hypothesis of no linear correlation (r = 0). Effective population size (N e ) estimates. The pruned set of 9,212 SNPs was used for the calculation of effective population size (N e ). To estimate N e , plink-formatted data was first converted to GENEPOP format using PGDSpider 2.1.1.5 87 . Then the LD method in NeEstimator2x 88 was used to calculate N e using the converted GENEPOP data as input. N e was calculated for global thoroughbreds and individual regions (ANZ, EUR, NAM and SAF). Year of birth for the horses used for N e calculation were restricted to 2013-2017.
Composite selection signals (CSS). The composite selection signals (CSS) method following the procedure of Randhawa et al. 24,25 was used to identify signatures of selection in elite Thoroughbred horses (n = 110) using the Thoroughbred and Founders data set summarised in table two. A mixed set of non-Thoroughbred horses representing putative founder populations (MOR, n = 23; CMP, n = 17; ARR, n = 24; AKTK, n = 20) as the comparator population. The Elite Thoroughbred population was assigned as the population under selection and the non-Thoroughbreds were assigned as the reference population.
The CSS approach was developed to investigate genomic signatures of selection and has been successful at localizing genes for monogenic and polygenic traits under selection in livestock 24,25,89 . The CSS uses fractional ranks of constituent tests and does not incorporate the statistics with P values, allowing a combination of the evidence of historical selection from different selection tests. For the present study, the CSS combined the fixation index (F ST ), the change in selected allele frequency (ΔSAF) and the cross-population extended haplotype homozygosity (XP-EHH) tests into one composite statistic for each SNP. F ST statistics were computed as the differentiation index between the population/s of interest (i.e. selected) and the contrasting/reference population/s (i.e. non-selected). XP-EHH and ΔSAF statistics were computed for the selected population(s) against the reference population. The CSS were computed as follows: (1) For each constituent method, test statistics were ranked (1,…, n) genome-wide on n SNPs.
(3) Fractional ranks were converted to z-values as z = Φ−1(r´) where Φ−1(⋅) is the inverse normal cumulative distribution function (CDF). (4) Mean z scores were calculated by averaging z-values across all constituent tests at each SNP position and P-values were directly obtained from the distribution of means from a normal N (0, m -1 ) distribution where m is the number of constituent test statistics. (5) Logarithmic (-log 10 of P-values) of the mean z-values were declared as CSS and were plotted against the genomic positions to identify the significant selection signals. (6) To reduce spurious signals, the individual test statistics were averaged (smoothed) over SNPs across chromosomes within 1 Mb sliding windows.