Some maternal lineages of domestic horses may have origins in East Asia revealed with further evidence of mitochondrial genomes and HVR-1 sequences

Samples and DNA extraction

Blood samples of 31 horses were collected from 14 Chinese domestic breeds across northern, northwestern and southwestern China (Table S1). The total genomic DNA was extracted using the phenol/chloroform extraction method. The total genomic DNA samples were further diluted to 50 ng/µl for the polymerase chain reaction (PCR) and DNA sequencing. All experimental procedures and protocols for sample collections were approved by the State Key Laboratory for Agro-biotechnology of CAU (Permit Number: XK257).

PCR amplification and mitochondrial DNA sequences

To amplify the whole mitochondrial DNA sequences of the Chinese horses, primers were designed according to the published sequence X79547 (Xu & Arnason, 1994) (Table S2). All of PCR amplifications were conducted in a 20 µl volume containing 10 µl MasterMix (TIANGEN, Beijing, China), 8 µl double distilled water (ddH₂0), 0.5 µl each primer (20 pmol/µl), and approximately 50 ng DNA. The PCRs were carried out as follows: 10 min of initial denaturation at 95 °C, followed by 30 cycles of 30 s at 95 °C, 60 s at annealing temperature (Table S2), and 30 s at 72 °C, with a final extension of 10 min at 72 °C. PCR products were directly sequenced using ABI 3730xl Genetic analyzer (Applied Biosystems, Foster City, CA, USA) by BGI-Tech Company (Beijing, China). In cases where samples were not amplified in the first attempt, re-amplification was performed. Sequences were deposited in the GenBank and accession numbers were listed in Table S1 .

We evaluated 31 de novo mitochondrial genomes with 317 others from the GenBank plus 16 sequences belonging to 16 Przewalskii horses (Equus przewalskii). The related information of the whole mtDNA sequences, including the geographical distribution and breed/type (where applicable), was provided in Table S3.

The HVR-1 sequences (5,180 sequences) in our study include: (1) 31 HVR-1 sequences isolated from the Chinese indigenous horses sequenced in this study (Table S1); (2) 317 HVR-1 sequences truncated from the whole mtDNA sequences obtained from Genbank (Table S4); and (3) 4,832 sequences retrieved from GenBank (Table S5).

Sequencing data analysis

By using MitoToolPy (http://www.mitotool.org/mp.html) (Peng & Fan, 2015), we scored the variants for each sequence (whole mtDNA genome sequences) relative to the reference sequence JN398377 (Achilli & Olivieri, 2012). We assigned the mtDNA sequences into specific haplogroups according to DomeTree (http://www.dometree.org/trees/horse.htm) (Peng & Fan, 2015). We followed the caveats for quality control in mtDNA dataset of domestic animals (Shi & Fan, 2014).

Phylogenetic reconstruction

To depict the phylogeographic patterns for whole mtDNA sequence, we assigned each of the 348 whole mtDNA sequences into haplogroups in context of updated haplogroup tree according to DomeTree (http://www.dometree.org/trees/horse.htm) (Peng & Fan, 2015). This tree encompasses 348 sequences and was rooted with donkey (Equus asinus) sequence.

NETWORK 5.0 (http://www.fluxus-engineering.com) was used for calculation of the median joining network for HVR-1 sequences. Calculations were performed under the following conditions: (1) default setting of Epsilon (0) was chosen; (2) three mutational hotspot sites 15585,15597 and 15650 were set as zero, and two mutational hotspot sites 15659,15737 were downweighted (weight :5) (Jansen & Forster, 2002; Cieslak & Pruvost, 2010); (3) the “>1 frequency” setting in MJ calculation was activated and MP Calculation option was conducted to remove excessive links and median vectors; and (4) other parameter settings were maintained as defaults.

Haplotype partitioning was quantified by DnaSP v5.0 software (Librado & Rozas, 2009). To verify the results of Network, we also assigned each of the haplotypes of HVR-1 sequences into haplogroups in context of updated mitogenome haplogroup tree according to DomeTree (http://www.dometree.org/trees/horse.htm) (Peng & Fan, 2015).

Molecular dating

Using 348 mtDNA genomes, we inferred the time of the most recent common ancestor of the haplogroups as employed in BEAST v1.8 (Drummond & Rambaut, 2007; Drummond & Suchard, 2012). We used GTR+I+G substitution model selected by jModelTest 2.1.4 (Darriba & Taboada, 2012; Guindon & Gascuel, 2003). The mutation rate of 7.39 × 10⁻² mutations/site/Mya (95% confidence interval 2.49 × 10⁻²–1.60 × 10⁻¹) (Lippold & Matzke, 2011) was used under both strict and relaxed clock models (Drummond & Rambaut, 2005). Each Markov chain Monte Carlo was run for 60,000,000 generations and the first 10% of each run was discarded as burn-in, and sampled every 1,000 generations. We confirmed convergence with Effective Sample Size values >300 in Tracer v1.6 (Drummond & Rambaut, 2007). We then used TreeAnnotator v1.6.1 implemented in the BEAST for the phylogenetic tree summaries and viewed the final tree with attributes such as posterior, node age, and height median in FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

Haplogroup tree and network

The analysis of whole mtDNA sequences revealed 16 clades (Aw-Rw) according to DomeTree (http://www.dometree.org/trees/horse.htm) (Peng & Fan, 2015) (Fig. 1). The system of nomenclature of those haplogroups was in accordance with Achilli & Olivieri (2012). The result showed that horses from the same breeds and/or regions were widely distributed within the clades. Further, the resulting matrilineal genealogy does not show evidence of strong concordance with breeds and/or regions. This therefore may support the hypothesis that horse may have multiple origins of matrilineal inheritance.

The alignment of the 5,180 HVR-1 sequences produced 397 haplotypes. We identified three strong mutational hotspots which have been previously described at nucleotide positions 15585, 15597, and 15650 (Vilà & Leonard, 2001; Jansen & Forster, 2002). These positions were not included in our phylogenetic analyses. In addition, two mutational hotspots 15659, 15737 were individually down-weighted (weight 5). Once these mutational hotspots were eliminated from the dataset, we divided the sequences into nine haplogroups from A to I based on the structure of the network (Fig. 2), of which haplogroups H and I were identified in our previous study (Yang & Zhu, 2017), and the four clusters B3, D5, D6, and G1 were newly observed in the present study.

Based on the result of Network we found a small proportion of HVR-1 haplogroups, which were out of the scope of A-I haplogroups. For completeness those we assigned each of the haplotypes of HVR-1 sequences into haplogroups in context of updated mitogenome haplogroup tree according to DomeTree (http://www.dometree.org/trees/horse.htm) (Peng & Fan, 2015), results were showed in Fig. 3 and Table S5.

Corresponding relationships between haplogroups of mtDNA HVR-1 and whole genomes

To establish the relationships between the haplotypes defined from HVR-1 and whole mtDNA sequences, we retrieved the HVR-1 sequences from the 348 whole mitochondrial genomes and conducted phylogenetic analysis (Table S4). A clear corresponding relationship was revealed when comparing the haplotypes of HVR-1 and mitochondrial genomes of all the samples in Table 1(except three Przewalskii horse and one ancient Yakut horse). A haplogroup of the genome corresponded to one or more of its HVR-1 sequence counterparts. The haplogroup A is the clade containing the largest number of haplotypes defined on the basis of HVR-1, but its sub-haplogroups clustered into 6 haplogroups of mitochondrial genomes. L_w contained three sub-haplogroups of HVR-1 and is the largest haplogroup containing most haplotypes of the mitochondrial genome.

Geographic distributions of horse mtDNA haplogroups

The analysis of geographic distribution of the mitochondrial genome haplogroups showed that horse populations in Europe or East Asia included all haplogroups defined from the mtDNA genome sequences. The lineage F_w comprised entirely of Przewalskii horses. The two haplogroups I_w and L_w displayed frequency peaks in Europe (14.08% and 37.32%, respectively) and a decline to the east (9.33% and 8.00% in the West Asia, and 6.45% and 12.90% in East Asia, respectively), especially for Lw, which contained the largest number of European horses (Table 2). However, an opposite distribution pattern was observed for haplogroups Aw, Hw, Jw, and Rw, which were harbored by more horses from East Asia than those from other regions. The proportions of horses from East Asia for the four haplogroups were 38%, 88%, 62%, and 54%, respectively.

Considering the large amount of HVR-1 sequence data of horse mtDNA currently available in the literatures or deposited in GenBank, we also analyzed the geographic distributions of the HVR-1 haplogroups (Table 3). We further confirmed our geographic distributions of the whole mtDNA haplogroups by analyzing the known mitochondrial HVR-1 sequences, using the corresponding relationships between the two sets of haplogroups (Table 1). The results showed that haplogroups D1, D2, and D3, which correspond to haplotype L_w of the whole mitochondrial sequences, displayed the highest frequencies in European horses, while I(corresponding to J_w), G1 (corresponding to R_w) showed frequency peaks in East Asia and a decline from Asia to Europe (Table 3). However, the patterns of geographic distribution of Aw and Hw were not supported by HVR-1 data. It should be noted that Cw contains H and majority of haplotypes in the HVR-1 haplogroup A6, which had a frequency peak in the West Asian regions. It is in accordance with geographic distribution of haplogroup A6, which is the most ancient haplogroup defined by HVR-1 sequences and contains the highest proportion of ancient horse sequences. However, haplogroup H has the highest frequency in the horse populations of East Asia.

Divergence times of haplogroups C_w, J_w, and R_w

The results showed that modern horses shared a common ancestor around 99,800 years ago. The phylogenetic tree for the estimation of the divergence time had 95% credibility intervals for the main branches. The resulting tree topology (Fig. 4) was very similar to the results obtained with the haplogroup tree of mtDNA genome sequences. Almost half of the haplotypes formed before the widely regarded time of the earliest horse domestication (approximately 5,800 years ago) (Outram & Stear, 2009).

From the phylogenetic analysis of both mitochondrial HVR-1 and genome sequences, haplogroups J_w, R_w, and some maternal lineage of Cw (H haplogroup defined by HVR-1) showed high frequency in East Asia. To further investigate their possible origins, we determined their divergence time. The age of haplogroup C_w was estimated to be 7.8 thousand years (KYA; 95% highest posterior density (HPD): 1,000–18,800 years), while the age of haplogroup J_w was 13.1 thousand years (KYA; 95% highest posterior density (HPD): 1,700–32,600 years) and haplogroup R_w has an age more than 6.7 thousand years (KYA; 95% highest posterior density (HPD): 900–16,600 years)(Table S6).