Genetic variability of Akhal-Teke horses bred in Italy

DNA extraction, genotyping and sequencing

The Italian AKH Association’s reference laboratory for genetic testing for parentage control at the Department of Veterinary Medicine (University of Milan) maintains the bio-bank of the biological samples of all AKH horses bred in Italy since they were first imported (107).

All applicable international, national and/or institutional guidelines for the care and use of animals were followed. The bio-banked samples were approved by the Ethics Committee of Università degli Studi di Milano Prot. OPBA_56_2016.

Ninety-five of the 107 blood samples (males and females) were genotyped using 17 microsatellite markers of StockMarks for the Horses^® Genotyping Kit (utilized for parentage testing), including 12 from the official ISAG panel (http://www.isag.us/Docs/EquineGenParentage2016.pdf). All PCR products were run on ABI PRISM^® 310 Genetic Analyzer, and fragments were analyzed with the Genescan^® and Genotyper^® software (http://tools.thermofisher.com/content/sfs/brochures/cms_040300.pdf).

According to the available pedigrees (first generation/parents), the horses with no common mare and stallion (n = 22) were sequenced, and so were 11 additional Arabian horses reared in Italy and bio-banked at the University of Milan. The 379 bp fragment of the mtDNA D-loop region (nucleotide position: 15,382–15,778) was amplified and sequenced by classic Sanger methods, using previously published primer pairs (Cozzi et al., 2004). The PCR products were purified utilizing SEPHADEX^® G-50 (SIGMA) and sequenced for both strands by using the Big Dye^® Terminator v.3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) according to the manufacturer’s guidelines. The PCR products were processed on an ABI PRISM^®310 Genetic Analyzer, and the nucleotide sequence was determined using the ABI PRISM^® DNA Sequencing Analysis Software.

Statistical analysis on microsatellite data

Allelic frequencies, genetic equilibrium according to Hardy–Weinberg (HW) (P-val), linkage disequilibrium (LD) and inbreeding coefficient (F_IS) were estimated using the GENEPOP software (Raymond & Rousset, 1995). LD between a pair of loci was also tested for genotypic data using the likelihood-ratio test of Arlequin v.3.5.2.2 software (Excoffier & Lischer, 2010). The number of alleles (NA), the effective number of alleles (NE) and the observed (H_o) and expected (H_e) heterozygosity per locus and on the entire population were calculated with POPGENE v.1.32 (Yeh, Yang & Boyle, 1997). Allelic richness (AR) per locus and in the population was also used to estimate genetic population diversity. AR was standardized for variations in sample size, and was calculated using FSTAT 2.9.3 (Goudet, 2002), based on 79 diploid individuals with no missing genotypes.

Statistical analysis on sequence data

Sequences were checked for the presence of ambiguous bases by using the software Chromas v.2.5.1 (http://www.technelysium.com.au/).

In order to define the haplotypes (H-) all obtained sequences were aligned with the complete horse mtDNA sequences X79547 (Xu & Arnason, 1994) and JN398377 (Achilli et al., 2012), widely used as reference sequences, using the ClustalW algorithm implemented in BioEdit (Hall, 1999).

The number of polymorphic sites, the parsimony informative and the singleton site, the number of haplotypes, the haplotype diversity (hd), the nucleotide diversity (π) and the average number of nucleotide differences (k) were estimated using the DnaSP5 v.5.10.01 software (Librado & Rozas, 2009). Haplotype pairwise genetic differentiation was calculated using Arlequin v.3.5.2.2 (Excoffier & Lischer, 2010).

In order to define the genetic relationships among our mtDNA sequences of AKH and those available in Genbank for other AKH lines and for other populations originated from the same Middle Eastern clade, a total of 104 sequences (two reference mtDNA, two Przewalski, five Arabian horse, 15 Barb horse and 13 Caspian pony, 56 Akhal-Teke and the 11 Arabian horses bred in Italy) were considered. All the sequences were trimmed at 254 bp to maximize the number of sequences included. Haplogroup definition was performed according to Achilli et al. (2012). The mtDNA sequences from GenBank database used for haplogroup definition and relationship investigation were reported in Table S1.

Three different approaches were used to investigate the haplotype and haplogroup relationship:

1. the Median-Joining Networks (MJN) analysis performed by Network v.5.0 software (http://www.fluxus-engineering.com); the default setting of Epsilon (0) was applied and the mutational hot spots present in the 254bp fragment were down-weighted as suggested by Jansen et al. (2002) and Cieslak et al. (2010);

2. the Neighbor-joining (NJ) tree (Saitou & Nei, 1987) based on the Kimura-2-parameter model distances (Kimura, 1980) was calculated among the haplotypes with MEGA v.7.0.14 software (Kumar, Stecher & Tamura, 2016); the bootstrap analysis, running 1,000 bootstrap replicates, was applied to check the robustness of the resulting dendrogram and a cut-off bootstrap at a significant value of 50 was applied;

3. the Principal Component Analysis (PCA) with a minimum spanning tree based on haplotype sequences was calculated using the PAST v.2.7 software (Hammer, Harper & Ryan, 2001).

The pairwise population genetic differentiation (F_ST) was calculated by Arlequin v.3.5.2.2 (Excoffier & Lischer, 2010) and used to evaluate genetic differentiation among Italian AKH, Italian Arabian, and the 93 horse GenBank sequences.

All the horse mtDNA sequences obtained in this study were deposited (GenBank ID: MF580226–MF580244).

Microsatellite markers

The results of microsatellite analyses are reported in Table 1. A total of 112 allelic variants were found on 17 microsatellites, ranging from four (HTG6 and HTG7) to 10 alleles (ASB2). The microsatellite LEX03 (nine alleles), mapping on X chromosome, was excluded from the statistical analyses.

The mean number of alleles (NA) calculated on 16 microsatellites was 6.4, whereas the effective mean number of alleles (NE) was 3.1. The analyses of the Hardy–Weinberg genetic equilibrium (P-val), calculated using HW formulae, showed three microsatellites (AHT4, HTG4 and CA425) that were not in genetic equilibrium (P < 0.05). This drift is probably caused by strongly unbalanced allelic frequencies in the small population originally imported to Italy from 1991 to 2000 (24 horses: 17 females and seven males).

According to the position (cM) of the microsatellites as reported by Swinburne et al. (2006), LD analyses were performed only for the markers located on the same chromosome: HMS6 and HTG7 (chr 4), HMS3 and HTG4 (chr 9) and ASB2, HTG6 and HMS1 (chr 15) (Table S2). LD between a pair of loci was tested in all samples together and by groups of individuals: those imported in Italy from 1991 to 2000 (Founders) and the ones born in Italy after 2000 (New population). Considering all samples together, all pairs of loci were in LD, except for HTG7 vs HMS6 (Table S2). The group of Founders was in equilibrium for all pair of loci; instead, for the New population the LD was observed for ASB2 vs HMS1 (GENEPOP and Arlequin) and ASB2 vs HTG6 (Arlequin) pairs of markers (Table S2). Hence, we calculated the NE based on the observed allelic frequencies that ranged from 1.7 (HTG4) to 5.3 (ASB2), thus revealing that some loci had a very unbalanced representation of allele variants. These included AHT4 and HTG4, two of the three loci that were not balanced, which present an NE value of 2.6 and 1.7, respectively.

Allelic richness (AR) is a parameter that is used to estimate diversity in the population. It ranges from 4 (HMS1) to 9.8 (ASB2), with a mean value of 6.4, comparable to the AR of 6.09 reported by Warmuth et al. (2011).

H_o and the H_e values ranged from 0.379 (HTG4) to 0.810 (ASB2), with a mean value of 0.642 and, from 0.410 (HTG4) to 0.816 (ASB2), with a mean value of 0.649, respectively The F_IS values per locus that ranged from −0.113 (HTG10) (outbreeding) to 0.174 (HTG6) (inbreeding); the overall F_IS value for the population was 0.016.

We estimated the diversity parameters on the AKH horses imported in Italy from 1991 to 2000. The values obtained were compared with those estimated in the AKH population born in Italy (after 2000). The two groups of horses are in HW equilibrium and have similar genetic values, except for the F_IS value that denotes, in the first nucleus, an apparent excess of heterozygosis Table S3.

Mitochondrial DNA sequences

Table 2 shows the eleven haplotypes found in AKH reared in Italy (H1 to H11). The 22 polymorphic sites identified were 19 parsimony informative and three singleton sites (15495, 15538 and 15709). All the polymorphic sites reported in Table 2 were transition sites. The maximum composite likelihood estimate of the nucleotide substitution pattern, identified for the 22 sequences, were 32.5% (A), 26.6% (T), 29.1% (C), and 11.8% (G), depending on the order of nucleotide composition in the vertebrate mitochondrial genome (A>C>T>G). The A+T content at the mtDNA D-loop region was 59.1% and exceeded G+C content (40.9%). Zhang et al. (2012) for Chinese Native horse breeds and Sziszkosz et al. (2016) for Hungarian Gidran horses report similar rich A+T content values in horse mtDNA, precisely 55.8% and 64.1% respectively.

Nucleotide differences compared to reference sequences ranged from 1 (H1) to 11 (H11). The most frequent haplotype found in our samples was H4 (18.2%), while the less frequent haplotypes were H1, H6, H8, H10 and H11 recorded in only one individual each (4.5%).

The values of variability parameter haplotype diversity (hd), nucleotide diversity (π) and average number of nucleotide differences (k) were 0.938, 0.021 and 6.448, respectively.

The pairwise distance matrix (Fig. 1) represents the number of nucleotide substitutions in each haplotype, where Hap_R represent the haplotype of the reference sequences. Inter-haplotypic differences ranged from 1 (Hap_R vs Hap_1 and Hap_10 vs Hap_11) to 13 (Hap_4 vs Hap_10 and Hap_4 vs Hap_11).

A total of 51 haplotypes with 36 variable sites (seven singleton and 29 parsimony informative) were identified by the alignment of the 126 mtDNA D-loop sequences, showing an hd value of 0.971. The detailed sequence variations of haplotypes and the identified haplogroups were reported in Table S4.

Four haplotypes (H16, H18, H21 and H22) were found exclusively in Italian AKH and were not reported in GenBank for the same breed. Additionally, our AKH sequences were grouped in five haplogroups A, C, G, L and M, as defined by Achilli et al. (2012).

The MJN (Fig. 2A) showed the relationships among the 51 haplotypes occurring in the 126 mtDNA sequences, and discretized them in the eight haplogroups. Relationships among sequences were also defined in the NJ tree (Fig. 2B) and in the PCA analyses (Fig. 2C) that, nevertheless, are unable to dissect the variation among haplogroups as well as the MJN. Finally, the matrix of pairwise population genetic differentiation (F_ST) showed a remarkable similarity between our AKH and other AKH populations, as expected (Table S4).

The NJ tree based on Kimura-2-parameter model distances and performed by applying a cut-off value of 50 for bootstrap value significance clearly separated haplogroups L from the others (62% of bootstrap value) (Fig. 2B). Haplogroups A and Q are divided into a different branch of the dendrogram (Fig. 2B). Particularly, haplotypes H2, which belong to haplogroup A, group together with haplotypes belonging to haplogroups C in another branch of the tree. The H2 located in the same operational taxonomic unit (OTU) of H3 and H50 differs from these by one mutation 15,650 (H2 vs H3) and by two mutations 15,650 and 15,728 (H2 vs H50), respectively. The haplotype H42, formally belonging to haplogroups G in the PCA but separated with haplogroup L in NJ tree, will be considered as intermediate among haplogroups G and L as well defined in MNJ (Figs. 2A, 2B and 2C).

The PCA analysis showing a minimum spanning tree based on sequence haplotypes, is displayed in Fig. 2C. The graphic placed in Fig. 2C reports the percentage of variance explained by PCA components. The first component (PC1) that explains most of the variability (21.2%) clearly separated haplogroups G+F and L, whereas the second component (PC2) (17.9%) defines haplogroup Q, confirming the MJN and NJ tree evidence. Haplogroups A, C, M and P appear mixed and are not mutually well separated both on PC1 and on PC2. The PCA analysis confirms the central role of H2.