Identification of an ancestral haplotype in the mitochondrial phylogeny of the ovine haplogroup B

Sample collection

All the animal procedures were performed in strict accordance with the guidelines of the Ethics Committee of Sassari University, Italy, which also approved this study. The biological samples were obtained from a female mouflon which died for natural causes, belonging to a feral population living in the Montes forest (40°12′19.08″–N 9°21′16.02″E), located within the Gennargentu national park in the mountainous Central-Eastern part of Sardinia. This population remained isolated from contacts with domestic species up to the present day and is considered, from a historical point of view, one of the oldest Sardinian mouflon colonies.

DNA extraction, amplification and sequencing

Genomic DNA was extracted from muscle using the GenElute Mammalian genomic DNA miniprep kit (Sigma-Aldrich) according to the manufacturer’s protocol. Sample quality and DNA concentration were determined via spectrophotometry using a ND-8000 (NanoDrop Technologies, Thermo Fisher Scientific Inc., Wilmington, DE, USA). Total yield of DNA extracted was ∼21 µg (DNA conc. 125–150 ng/µl with an elution volume of 150 µl).

The 21 pairs of primers reported in Sanna et al. (Sanna et al., 2015) were used to amplify overlapping fragments spanning the whole mitogenome of the Sardinian mouflon.

PCR products were purified by electrophoresis on a 1.2% agarose gel using the Montage™ DNA Gel extraction kit (Millipore) and then sequenced with each of the PCR primers on an ABI 3130 sequencer. Sequences obtained were aligned, edited and assembled using Clustal X2 (Larkin et al., 2007) and the software CAP3 Assembly Program (Huang & Madan, 1999) to generate the complete mitogenome sequence. The protein-coding genes were determined by open reading frame finder implemented at the NCBI website with the vertebrate mitochondrial genetic code, and then finally confirmed by sequence comparisons with the reported O. aries mitogenome (GB# NC_001941).

The physical map of the Sardinian mouflon mitogenome was generated by means of OGDraw 1.2 (Lohse et al., 2013).

Homologous sequences assembly, dataset composition and phylogenetic analyses

Three data set including 69 D-loop (Table S1), 90 CytB (Table S2) and 29 complete mitogenome (Table S3) sequences, respectively, were assembled for phylogenetic analyses.

For the complete mitogenome, the comparisons were performed on both the entire mtDNA sequences and 28H marker encompassing the concatenated sequences of the 12 protein-coding genes, 14 transfer RNA genes and two ribosomal RNA genes, all located on the H-strand, as previously described (Sanna et al., 2015). We produced a joined dataset using the complete mitogenome sequence of a Sardinian mouflon from the Montes area obtained in this study, plus 27 sequences collected from GenBank, and the mitogenome sequence of an Asian mouflon (KF938360), recently characterized (Lv et al., 2015) (Table S3). Ingroup taxa included 17 Ovis species, the Moschus moschiferus was selected to root the phylogenetic trees while the remaining 11 outgroup sequences were used as calibration points for the divergence time estimations.

Sequences were aligned using the Clustal X2 software (Larkin et al., 2007). DNA sequence variation parameters were estimated using the software package DnaSP 5.10 (Librado & Rozas, 2009). The best-fit model for each dataset was selected using MEGA 7.0.14 (Kumar, Stecher & Tamura, 2016). The general time-reversible model with discrete gamma distribution and invariant sites (GTR+G+I) had the highest score for both the D-loop and the complete mitogenome sequences, while GTR+G was the best for the CytB data set.

Phylogenetic relationships among individuals were investigated by means of Bayesian Inference (BI) and Maximum Likelihood (ML) analyses. MrBayes 3.2.4 (Ronquist et al., 2012) was used for BI and MEGA 7.0.14 (Kumar, Stecher & Tamura, 2016) for ML tree reconstruction. Analyses for BI tree reconstruction were carried out for five million generations and the trees were sampled every 10 generations. The reliability of each branch in ML tree was estimated by bootstrapping (10,000 replications)

The TCS network analysis was performed using the software Pop ART 1.7 (Leigh & Bryant, 2015), to identify possible disconnections between groups of individuals, and infer the genetic relationships among the haplotypes.

Divergence time estimation was inferred by the complete mitogenome sequences analysis by assuming seven different calibration points (CPs) based on fossil records providing ages for nodes inside Bovidae (Bibi, 2013) as previously described (Sanna et al., 2015).

Mitogenome characterization

We sequenced the complete mtDNA sequence of a Sardinian mouflon individual and deposited it in GenBank (MG489885). The obtained 16,618 bp long molecule showed a percentage composition in bases of the L-strand of 33.7% A, 25.8% C, 13.1% G, and 27.4% T. Overall, we detected 13 protein coding genes, 12S and 16S ribosomal RNA, 22 transfer RNA genes and a D-loop region, according to the molecule gene content of the other Ovis species previously reported (Hiendleder et al., 1998; Lv et al., 2015; Meadows, Hiendleder & Kijas, 2011; Sanna et al., 2015). The structural organization of Sardinian mouflon mitogenome and the location of the various features are displayed in Fig. S1 and summarized in Table 1.

D-loop and CytB sequences analysis

A preliminary phylogenetic analysis was performed using the mtDNA D-loop region and the CytB gene as molecular markers. The data set of D-loop region consisted of 69 sequences from domestic sheep and 17 from wild mouflons (Table S1). The urial sequence was selected as outgroup to root the BI tree (Fig. 1).

Five domestic sheep HPGs were detected, accordingly with those described by Meadows, Hiendleder & Kijas (2011) on the entire mitogenome sequence analysis. Three main evolutionary lineages within the sheep/mouflon early radiation gave birth to HPG D, the clade including HPGs A and B, and another consisted of HPGs C and E grouped together with haplotype (Hpt X), respectively. The sequences from wild species were not uniformly distributed throughout the HPGs. The European mouflon sequences (OAMs in Table S1) all belong to HPG B where even a sequence from an Asian mouflon (OG) falls.

The OAM sequences were further grouped in three distinct clades. The OAM_3 sequence was not included within any group. The ML analysis (Fig. S2) provided an identical but slightly lower supported tree topology.

To further investigate the relationships among wild and domestic species within the HPG B, a TCS network analysis was carried out (Fig. 2).

By contrast to the other European and Asian mouflon sequences the position of the OAM_26 and OAM_5 could support an early split from the urial lineages before the advent of the first domestic sheep species originating from a central common ancestor (kk2 nucleus in Fig. 2).

The Cyt B sequences analysis allowed to significantly increase the number of wild species included in the investigation (Table S2). Indeed, the O. orientalis sequences filled the 34.4% of the entire dataset, 6.7% were from hybrid O. orientalis/O. vignei individuals, 7.8% from O. vignei and the remaining 50% belonged to O. aries. The OAM_26 sequence was the only one to be selected from the Sardinian group because of its crucial position, together with the OAM_5, in the D-loop analysis (Fig. 2). The BI tree analysis (Fig. 3) evidenced the raise of four mtDNA lineages within the O. aries/O. orientalis clade: three already detected in the D-loop sequences analysis and a further one originating the OOG4-OOG5 sequences group.

The results were similar to those obtained in D-loop analysis also including an increased number of wild species: all the European mouflon sequences are grouped within the HPG B, where no Asian mouflon sequence falls, except for the OG one, which given the matrilineal transmission of the mtDNA, could be sampled from a hybrid individual. Moreover, the OAM_26 sequence was found to be not related to sequences from German mouflons, as indicated in the TCS network (Fig. 4).

The ML analysis (Fig. S3) provided an identical but slightly lower supported tree topology.

Complete mitogenome phylogeny and molecular dating

Phylogeny on whole mitogenome was inferred using a Sardinian mouflon sequence from this study (MG489885), the dataset from Sanna et al. (2015), and an Asian mouflon sequence collected from GenBank (KF938360) (Table S3). The BI tree analysis performed on both the whole mitogenome and the 28H dataset supported the same relationships among species as confirmed by the identical tree topologies obtained (Fig. 5).

The ML trees (Fig. S4) obtained from the two datasets were comparable with BI results, albeit showing slightly lower supported tree topology. We identified two main clusters in the sheep/mouflon radiation (Prob ≥ 0.99) (Fig. 5).

Cluster II was characterized by two main groups: a first one including HPG-D, and a second one grouping HPGs A and B, and the Sardinian, mainland Europe and Asian mouflon. The Sardinian mouflon haplotype split before the rise of HPG-B sheep, Asian and mainland Europe mouflon, the latter two sharing a most recent common ancestor (MRCA).

The TCS network analysis (Fig. 6) of the whole Ovis mitogenome sequences retrieved three main disconnected clusters within the Ovis genus (Figs. 6A–6C).

Cluster A was exclusive of the wild sheep species included in the network analyses (argali, bighorn and urial). Cluster B encompassed HPGs D, A and B domestic sheep, along with Sardinian, mainland Europe and Asian mouflon sequences. Cluster C grouped the Cypriot mouflon and HPGs C and E domestic sheep. Within cluster B, the mainland Europe and Asian mouflon sequences constituted a sub-clade with HPG-B domestic sheep, while the Sardinian mouflon held an external position, showing a longer evolutionary process as confirmed by the high number of nucleotide substitution in its own lineage. A similar result was observed in the relationship among species within Cluster C, where a sub-clade encompassing HPGs C and E was detected, and the Cypriot mouflon held an external position.

The divergence time estimations dated back the early domestic sheep radiation to ∼410,000 YA, while the speciation events which gave rise to Cluster I and II dated 190,000 and 330,000 YA, respectively. Within Cluster II, the Sardinian mouflon haplotype split about 110,000 YA from the sub-clade including HPG-B sheep and mainland Europe mouflon sequence, which shared a common ancestor 8,000 YA. According to the results provided from the previous analyses carried out on the D-loop e CytB sequences, the KF938360 sequence was grouped within this same sub-clade, something that confirms the hypothesis that it may have been collected from a hybrid individual.

Genetic variability estimation based on mitogenome sequence analysis

Among the 17 complete mitogenomes belonging to five species within the Ovis genus, 17 haplotypes defined by 1,387 polymorphic sites (S) were found. Overall, high values of total mean haplotype and nucleotide diversity were obtained (h = 1 and π = 0.016, respectively). A lower value of nucleotide diversity (π = 0.008) was found among the O. aries individuals.

To quantify the divergence among HPGs, pairwise genetic distances between groups were calculated under the K2P model. The variation rate among sites was modelled with a gamma distribution (shape parameter = 0.64). Among domestic sheep HPGs the average value of genetic differentiation was high at 0.0086 (± 0.0006); the lowest level was observed between HPGs C and E (0.0036 ± 0.0004), closely followed by HPGs A and B (0.006 ± 0.0055). The distance between wild species and domestic sheep HPGs ranged from 0.055 ± 0.0009 (Bighorn vs domestic HPGs) to 0.0072 ± 0.004 (European mouflon vs. domestic HPGs). Among mouflon species, the Sardinian mouflon showed the largest genetic distance from domestic sheep HPGs (0.0136 ± 0.004), followed by Cypriot (0.0087 ± 0.0028), Asian (0.00723  ± 0.004) and mainland Europe mouflon sequences (0.00717 ± 0.004).

The origin of the ovine mitochondrial HPG-B

HPG-B is the most common mitochondrial haplogroup among extant domestic sheep, including the vast majority of cosmopolitan sheep breeds (Hiendleder et al., 1998; Meadows, Hiendleder & Kijas, 2011). The molecular dating we performed shows the Sardinian mouflon haplotype (HPT-B1) stemming about 30,000 years earlier than the HPG-B clade (∼80,000 YA), which includes domestic sheep and mainland Europe mouflons (Fig. 2). We suggest the presence of multiple mitochondrial ovine HPTs populating the Fertile Crescent around the end of the last glacial period, including those identified in the Cypriot, Sardinian (HPT-B1), mainland Europe (HPT-B2), and Asian mouflon. Our results suggest HPT-B1 to be a putative ancestral haplotype of the HPG-B lineage. Among the extant wild species we analysed, HPT-B1was found exclusively in the Sardinian mouflon from the Montes area. However, analyses on an extended dataset comprising more mouflon individuals and populations (e.g., Corsican mouflon) are necessary to confirm this result.

Mouflon and HPG-B sheep colonization of mainland Europe

Two hypotheses for the predominance of HPG-B in European sheep have been proposed. The most recent and accredited one, which is strongly supported by zooarchaeological and molecular evidences, suggests the predominance of HPG-B in European sheep as the result of the founder effect driven by a subset of founder animals belonging to the HPG-B (Tapio et al., 2006). The second one invokes an independent European sheep domestication event (Ryder, 1981) but there is no zooarchaeological findings to supports this hypothesis.

Our results show that the Sardinian mouflon haplotype (HPT-B1) split 110,000 YA from the mainland Europe mouflon (HPT-B2) and HPG-B sheep (HTP-B3) lineages, the latter sharing a MRCA 80,000 YA (Figs. 2 and 3), confirming the first hypothesis. Likely, the first settlers happened to choose HPT-B1, HPT-B2, and HPT-B3 from a wider pool of haplotypes, and brought them along when they moved from the Near East to Western Europe after the last glaciations (14 to 12,000 YA) (Fu et al., 2016), eventually colonizing Europe, including Corsica and Sardinia (Chessa et al., 2009; Clutton-Brock, 1999). A fascinating hypothesis is that the mutation that modified the fleece from hirsute to woolly—which putatively occurred on HPT-B3 individuals—was key in determining the value of mouflon to their early farmers (Demars et al., 2017). Since the phenotypic traits related to the hair type are not related to mtDNA sequence modifications, the exclusivity of the HPT-B3 in the HPG-B domestic sheep would be the result of co-selection for non-related traits. Since then this trait underwent artificial selection and the economic interest for the mouflon-like ovines quickly decreased, leading the mouflon to survive as feral populations. Alternatively, only the Sardinian and possibly the mainland Europe mouflon haplotypes (HPT-B1 and B2) were included by chance within the first group of individuals moved to Europe from the Neolithic people, as part of what is commonly described as the “first wave of sheep domestication” (Chessa et al., 2009). Successively, individuals harbouring the HPT-B3—and selected for new and favourable secondary traits—arrived to Europe (second wave of domestication) and replaced the resident ovine populations (Chessa et al., 2009; Hiendleder et al., 1998; Poplin, 1979).