# Evolutionary history of Mexican domesticated and wild Meleagris gallopavo

Gabriela Padilla-Jacobo, Horacio Cano-Camacho, Rigoberto López-Zavala,

María E. Cornejo-Pérez, María G. Zavala-Páramo

► To cite this version:

Gabriela Padilla-Jacobo, Horacio Cano-Camacho, Rigoberto López-Zavala, María E. Cornejo-Pérez, María G. Zavala-Páramo. Evolutionary history of Mexican domesticated and wild Meleagris gallopavo. Genetics Selection Evolution, 2018, 50 (1), pp.19. 10.1186/s12711-018-0388-8. doi:10.1186/s12711-018-0388-6

HAL Id: hal-02405226

https://hal.science/hal-02405226

Submitted on 11 Dec 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Evolutionary history of Mexican domesticated and wild *Meleagris gallopavo*

Gabriela Padilla-Jacobo<fnref n="1" />, Horacio Cano-Camacho<fnref n="1" />, Rigoberto López-Zavala², María E. Cornejo-Pérez<fnref n="1" /> and María G. Zavala-Páramo<fnref n="1" />*

**Abstract**

**Background:** The distribution of the wild turkey (*Meleagris gallopavo*) extends from Mexico to southeastern Canada and to the eastern and southern regions of the USA. Six subspecies have been described based on morphological characteristics and/or geographical variations in wild and domesticated populations. In this paper, based on DNA sequence data from the mitochondrial D-loop, we investigated the genetic diversity and structure, genealogical relationships, divergence time and demographic history of *M. gallopavo* populations including domesticated individuals.

**Results:** Analyses of 612 wild and domesticated turkey mitochondrial D-loop sequences, including 187 that were collected for this study and 425 from databases, revealed 64 haplotypes with few mutations, some of which are shared between domesticated and wild turkeys. We found a high level of haplotype and nucleotide diversity, which suggests that the total population of this species is large and stable with an old evolutionary history. The results of genetic differentiation, haplotype network, and genealogical relationships analyses revealed three main genetic groups within the species: mexicana as a population relict (C1), merriami (C2), and mexicana/intermedia/silvestris/osceola (C3). Haplotypes detected in domesticated turkeys belong to group C3. Estimates of divergence times agree with range expansion and diversification events of the relict population of *M. gallopavo* in northwestern Mexico during the Pliocene-Pleistocene and Pleistocene-Holocene boundaries. Demographic reconstruction showed that an expansion of the population occurred 110,000 to 130,000 years ago (Kya), followed by a stable period 100 Kya and finally a decline ~ 10 Kya (Pleistocene-Holocene boundary). In Mexico, the Trans-Mexican Volcanic Belt may be responsible for the range expansion of the C3 group. Two haplotypes with different divergence times, MGMDgoB/MICH1 and MICH2, are dominant in domesticated and commercial turkeys.

**Conclusions:** During the Pleistocene, a large and stable population of *M. gallopavo* covered a wide geographic distribution from the north to the center of America (USA and Mexico). The mexicana, merriami, and mexicana/intermedia/silvestris/osceola genetic groups originated after divergence and range expansion from northwestern Mexico during the Pliocene-Pleistocene and Pleistocene-Holocene boundaries. Old and new maternal lines of the mexicana/intermedia/silvestris/osceola genetic group were distributed within the Trans-Mexican Volcanic Belt where individuals were captured for domestication. Two haplotypes are the main founder maternal lines of domesticated turkeys.

**Background**

*Meleagris gallopavo* is an original neartic species with a distribution that extends from Mexico to southeastern

Canada and to the eastern and southern regions of the USA [1, 2]. Six subspecies of *M. gallopavo* have been described based on their geographical distribution and morphological characteristics such as size, coloration or iridescence of plumage, color of the legs, and color of the tip and base of the feathers, i.e. (1) *M. g. gallopavo* (domesticated) described by Linnaeus in 1758, (2) *M. g. silvestris* (Silvestre) described by Vieillot in 1817, (3) *M. g. mexicana* (Gould) described by Gould in 1856, (4) *M.*

<footnote>*Correspondence: gzavpar@hotmail.com</footnote>

<footnote n="1">¹Centro Multidisciplinario de Estudios en Biotecnología, FMVZ,
Universidad Michoacana de San Nicolás de Hidalgo, Km. 9.5 Carretera
Morelia-Zinapécuaro, Posta Veterinaria, C.P. 58000 Tarímbaro, Michoacán,
Mexico</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
g. intermedia (Río Grande) described by Sennett in 1879, (5) M. g. osceola (Florida) described by Scott in 1890, and (6) M. g. merriami (Merriam) described by Nelson in 1900 (for more details see [2–5]). *Meleagris gallopavo* is the one and only important domesticated animal species of North American origin [6]. Molecular studies based on mtDNA have suggested that the domesticated turkey is representative of the extinct wild subspecies M. g. gallopavo [7, 8]. Knowledge from historical registries indicates that different prehispanic Mexican groups such as the Purepecha, Huicholes and other ethnic groups of domesticated turkeys were present between 200 and 700 BC [9]. Leopold [6] and Nelson [10] proposed that domestication occurred in the highlands of Michoacan, Mexico, and according to Schorger [3] domesticated turkey stocks were established by at least ca. 200 BC to 700 AD within the Tehuacan Valley (Puebla), with bones dated from approximately 700 AD being identified in Guatemala.

There are few molecular genetic studies on domesticated turkeys from rural Mexican communities. As far as we know, there is only one analysis that used microsatellite markers to analyze domesticated turkey populations from the five physiographic regions of Michoacan in Mexico and that revealed three genetically distinct groups [11]. Some studies focused on commercial or heritage turkeys, and for example, a microsatellite analysis showed that the commercial turkey is closer to the heritage Narragansett, Bourbon Red, Blue Slate turkeys than to the Spanish Black and Royal Palm turkeys [12]. A genomic study that included seven commercial lines, three samples of wild turkeys from Chihuahua in Mexico, and the heritage Beltsville Small White, Royal Palm and Narraganset varieties revealed that all commercial lines shared the same origin and that specific haplotypes may have been selected in the modern domesticated turkey [13]. Other studies have focused on the analysis of the diversity between subspecies and conservation of wild populations [5, 14, 15]. Finally, an analysis of samples from bones and coprolites from archaeological sites in the southwestern USA and from the six proposed M. gallopavo subspecies using mitochondrial markers proposed two sites of turkey domestication that each involve wild turkey populations i.e. (1) M. g. gallopavo in south-central Mexico and (2) M. g. intermedia/silvestris with a subsequent introduction of domesticated stocks in the southwestern USA [7].

For agriculturally important species such as chicken (Gallus gallus) [16–19], duck (Anas platyrhynchos) [20–23], cattle (Bos taurus) [24–26], and pigs [27, 28], phylogenetic and genealogical molecular analyses have helped to better understand the process of their domestication and to demonstrate the origin and inter- or intraspecific

relationships of these species. Moreover, the use of mtDNA sequences in phylogeography analyses has been extensively tested and offers a highly sensitive method to analyze evolutionary processes [29]. Currently, they are the most widely used markers for such studies in vertebrates.

In this study, our aim was to investigate the genetic diversity and structure, genealogical relationships, divergence times and the demographic history of M. gallopavo, by putting emphasis on domesticated individuals to reconstruct the evolutionary history of this species. For these analyses, we used sequences of the mtDNA D-loop from domesticated, commercial and wild turkeys.

## Methods

### Sample collection

Blood and tissue samples were collected from domesticated and wild turkey populations between 2001 and 2011. For each individual, blood samples of 0.1 to 0.2 mL were taken from the brachial vein or tissue fragments as large as a grain of rice were obtained and placed in 2 mL vials with 0.5 mL of storage and lysis buffer (100 mM Tris pH 8.0, 100 mM EDTA pH 8.0, 10 mM NaCl and 2% SDS), saturated with salt-DMSO for tissue samples [30]. Then, samples were stored at room temperature during transport (about 1 week) and subsequently stored at 4 °C until further analysis.

Samples were deposited in the Collection of Biological Samples of the Centro Multidisciplinario de Estudios en Biotecnología (CMEB) of the Universidad Michoacana de San Nicolas de Hidalgo. A total of 187 samples were available for analysis: 161 originated from domesticated populations of *M. gallopavo* collected in Mexican rural communities from different localities in Puebla and from the five physiographic regions in Michoacan (Bajio, Trans-Mexican Volcanic Belt, Balsas, Sierra, and Costa) [31] (see Fig. 1 and Table 1); four samples of domesticated turkeys were obtained from the Izabal department in the northeastern region of Guatemala; nine samples of the commercial line Bronze were collected on a farm in Ario de Rosales, Michoacan, Mexico (Fig. 1, Table 1); nine samples of *M. g. mexicana* individuals were obtained from a management unit for wildlife conservation in Canatlan, Durango in Mexico, a population that was reintroduced from Yecora, Sonora in Mexico; and four samples of *M. g. intermedia* were donated by hunters that held a permission to hunt in Villa de Casas, Tamaulipas in Mexico (Fig. 1, Table 1).

### DNA extraction and amplification

DNA was extracted from tissue and blood samples using the phenol-free method described by FitzSimmons [32]. The mtDNA D-loop sequence was amplified using
<figure><img src="image_3.png" /></figure>

the oligonucleotides NAU313 (5' GCCACCTGTGAA-GAAGCC 3') and NAU185 (5' ACGGCTTGAAAAGC-CATTGTTGT 3') [5]. PCR reactions were performed in a total volume of 25 µL as follows: PCR buffer 1× (20 mM Tris-HCl pH 8.4, 50 mM KCl, 1.5 mM MgCl₂), 200 mM of each dNTP, 10 pmol of each oligonucleotide, 1.5 U Platinum Taq polymerase (Invitrogen) and 50 ng of DNA. The reaction mixtures were placed in a thermocycler (Gene Amp 2700, Applied Biosystems) under the following amplification conditions: 95 °C for 5 min, followed by 30 cycles of 95 °C for 1 min, 60 °C for 1 min and

72 °C for 2 min, with a final extension at 72 °C for 8 min. DNA sequencing was performed using the dideoxy technique on both strands [33] using the commercial service Macrogen USA. In addition to these 187 sequences, we analyzed 425 sequences from wild and domesticated individuals that were obtained from the NCBI GenBank database (Table 2), which amounted to 612 sequences. Figure 1 shows the location of the individuals sampled for this study and the geographical origin of the sequences obtained from the NCBI GenBank database.
<table><caption>Table 1 Localities where samples of wild and domesticated turkeys from Mexico and Guatemala were collected</caption><thead><tr><th>Country</th><th>Sample types</th><th>State</th><th>Locality</th><th>Sample size</th></tr></thead><tbody><tr><td rowspan="15">Mexico</td><td>M. g. mexicana</td><td>Durango</td><td>Canatlan (UMA El Durangueño)</td><td>9</td></tr><tr><td>M. g. intermedia</td><td>Tamaulipas</td><td>Villa de Casas</td><td>4</td></tr><tr><td rowspan="12">Domesticated</td><td rowspan="12">Michoacán</td><td>Bajío Region (Penjamillo, Indaparapeo, Zacapu, La Piedad, Jimenez, Capula, Churintzio, Zamora)</td><td>33</td></tr><tr><td>Trans-Mexican Volcanic Belt Region (Zitacuaro, Cuitzeo, Benito Juarez, Tuxpan, Tacambaro, Ario de Rosales, Salvador Escalante)</td><td>31</td></tr><tr><td>Balsas Region (Churumuco, Huetamo, Tiquicheo, Tepalcatepec, Tuzantla, Nuevo Urecho, Buenavista, Apatzingan)</td><td>26</td></tr><tr><td>Sierra Region (Tumbiscatio, Arteaga, Aguililla)</td><td>28</td></tr><tr><td>Costa Region (Coahuayana, Aquila, Arteaga)</td><td>26</td></tr><tr><td rowspan="5">Puebla</td><td>Rafael Lara Grajales (Benito Juarez Ejido)</td><td>1</td></tr><tr><td>Totimehuacán (San Baltazar Azumiatla)</td><td>5</td></tr><tr><td>Soltepec (San Mateo)</td><td>7</td></tr><tr><td>Hueyotlipán (Sto. Tomas Hueyotlipán)</td><td>4</td></tr><tr><td>Ario de Rosales</td><td>9</td></tr><tr><td rowspan="2">Commercial line Bronze</td><td rowspan="2">Michoacán</td><td rowspan="2">Izabal</td><td rowspan="2">Izabal</td><td>4</td></tr><tr><td>4</td></tr><tr><td>Guatemala</td><td>Domesticated</td><td>Izabal</td><td>Izabal</td><td>187</td></tr><tr><td>Total</td><td></td><td></td><td></td><td></td></tr></tbody></table>

## Genetic diversity and differentiation

Sequence editing, alignment, and construction of data matrices were carried out with Sequencher v4.1 [34] and PhyDE [35]. The number of haplotypes (H), polymorphic sites (S), and nucleotide (π) and haplotype (Hd) diversity estimates for the domesticated and wild populations were calculated with DnaSP v5 [36]. Analysis of molecular variance (AMOVA) [37] was used to calculate genetic variation and genetic differentiation between populations by performing 10,000 permutations. In addition, computed pairwise comparisons of FST values with 1000 permutations were obtained with ARLEQUIN v3.1 [38].

## Genealogical relationships between haplotypes

To establish genealogical relationships between haplotypes and their frequencies, a haplotype network was constructed using the median-joining method [39] with the software NETWORK v4.6.0.0 [40] and setting default parameters. The relationships between haplotypes were also analyzed using phylogenetic inference. Matrices for these analyses included haplotypes that were identified in this study and haplotypes for each subspecies that are reported in the NCBI GenBank database (Table 2). The sister species *Meleagris ocellata* was included as outgroup (Table 2). Models of molecular evolution were evaluated with jModelTest v2.1.1 [41] and selected using the corrected Akaike Information Criterion (cAIC) [42]. The best model obtained using this criterion was Hasegawa, Kishino and Yano [43] + Invariant sites, i.e. HKY+I.

Reconstructions of genealogical relationships were generated using maximum likelihood (ML) and Bayesian inference (BI) frameworks with RAxML v8 [44] and MrBayes v3.2 [45], respectively. Branch support values were estimated by bootstrap analysis (BP) of 500 replicates and by calculating posterior probabilities (PP). MrBayes was run with the following parameters: four independent runs of four chains each (one cold chain and three hot chains) for 10 million generations, sampling one tree every 1000 generations. Trees and parameters were summarized after discarding 25% of the data as burn-in. The remaining trees were summarized as a majority consensus tree and visualized using FigTree v1.4.0 [46].

## Estimation of divergence times and rates of molecular evolution

The data matrix included two sequences of *M. ocellata* and one sequence of *Gallus gallus* (GenBank Access HQ022888.1). Divergence times were estimated using BEAST v1.7.4 [47]. An uncorrelated lognormal relaxed clock model was selected with the HKY+I model of evolution. Because of the nature of the data, a tree prior with a coalescent model assuming a constant population size was used [48]. One calibration point with a lognormal prior distribution (mean = 0.0, standard deviation (SD) = 1.0, offset = 2.6) and the oldest *M. gallopavo* fossil (2.6 Mya) that is registered in the PaleozoDB fossil database (82,258) [49] were used. Markov chain Monte Carlo (MCMC) analyses were run for 10 million generations,
<table><caption>Table 2 Nucleotide sequences of D-loop from *M. gallopavo* and codes used in the analyses</caption><thead><tr><th>Population/sub-species</th><th>GenBank accession number</th><th>Code in this study</th><th>Population/subspecies</th><th>GenBank accession number</th><th>Code in this study</th></tr></thead><tbody><tr><td rowspan="2">M. g. mexicana</td><td>AF486960</td><td>MGM60</td><td rowspan="2">M. gallopavo, archaeological samples</td><td>GQ303159</td><td>MGarq59</td></tr><tr><td>AY037888</td><td>MGM88</td><td>GQ303160</td><td>MGarq60</td></tr><tr><td rowspan="15">M. g. silvestris</td><td>AF486901</td><td>MGS01</td><td rowspan="5">M. g. gallopavo 1903</td><td>GQ303161</td><td>MGarq61</td></tr><tr><td>AF486902</td><td>MGS02</td><td>GQ303163</td><td>MGarq63</td></tr><tr><td>AF486905</td><td>MGS05</td><td>GQ303164</td><td>MGarq64</td></tr><tr><td>AF486911</td><td>MGS11</td><td>GQ303165</td><td>MGarq65</td></tr><tr><td>AF486914</td><td>MGS14</td><td rowspan="10">M. g. gallopavo</td><td>EF153719</td><td>MG19</td></tr><tr><td>AF486918</td><td>MGS18</td><td>AF172952</td><td>MG52</td></tr><tr><td>AF486921</td><td>MGS21</td><td>AJ297180</td><td>MG80</td></tr><tr><td>AF486927</td><td>MGS27</td><td>AF487103</td><td>MGI03</td></tr><tr><td>AF172947</td><td>MGS47</td><td>AF487117</td><td>MGI17</td></tr><tr><td>AF172953</td><td>MGS53</td><td>AF487058</td><td>MGI58</td></tr><tr><td>AF172954</td><td>MGS54</td><td>AF487059</td><td>MGI59</td></tr><tr><td>AF172957</td><td>MGS57</td><td>AF487060</td><td>MGI60</td></tr><tr><td>AF172958</td><td>MGS58</td><td>AF487062</td><td>MGI62</td></tr><tr><td>AF172960</td><td>MGS60</td><td>AF487063</td><td>MGI63</td></tr><tr><td>AF172961</td><td>MGS61</td><td>AF487065</td><td>MGI65</td></tr><tr><td>AF486875</td><td>MGS75</td><td>AF487067</td><td>MGI67</td></tr><tr><td>AF486876</td><td>MGS76</td><td>AF487071</td><td>MGI71</td></tr><tr><td>AF486877</td><td>MGS77</td><td>AF487072</td><td>MGI72</td></tr><tr><td>AF486885</td><td>MGS85</td><td>AF487077</td><td>MGI77</td></tr><tr><td>AF486887</td><td>MGS87</td><td>AF487085</td><td>MGI85</td></tr><tr><td>AF486889</td><td>MGS89</td><td>AF487094</td><td>MGI94</td></tr><tr><td rowspan="10">M. g. oceola</td><td>AF486895</td><td>MGS95</td><td rowspan="2">M. g. merriami</td><td>AF487006</td><td>MGMer06</td></tr><tr><td>AF486898</td><td>MGS98</td><td>AF487010</td><td>MGMer10</td></tr><tr><td>AF486931</td><td>MGO31</td><td>AF487023</td><td>MGMer23</td></tr><tr><td>AF486938</td><td>MGO38</td><td>AF487025</td><td>MGMer25</td></tr><tr><td>AF486944</td><td>MGO44</td><td>AF487039</td><td>MGMer39</td></tr><tr><td>AF486949</td><td>MGO49</td><td>AF487041</td><td>MGMer41</td></tr><tr><td>AF486951</td><td>MGO51</td><td>AF487042</td><td>MGMer42</td></tr><tr><td>AF486953</td><td>MGO53</td><td>AF172948</td><td>MGMer48</td></tr><tr><td>AF486956</td><td>MGO56</td><td>AF172964</td><td>MGMer64</td></tr><tr><td>AF486959</td><td>MGO59</td><td>AF486985</td><td>MGMer85</td></tr><tr><td rowspan="5">M. gallopavo, archaeological samples</td><td>GQ303154</td><td>MGarq54</td><td rowspan="5">M. ocellata</td><td>AF486986</td><td>MGMer86</td></tr><tr><td>GQ303155</td><td>MGarq55</td><td>AF486997</td><td>MGMer97</td></tr><tr><td>GQ303156</td><td>MGarq56</td><td>AF486999</td><td>MGMer99</td></tr><tr><td>GQ303157</td><td>MGarq57</td><td>AF487120</td><td>AF487120Aocellata</td></tr><tr><td>GQ303158</td><td>MGarq58</td><td>AF487121</td><td>AF487121Aocellata</td></tr></tbody></table>

sampling one tree every 1000 generations. The results were summarized using TreeAnnotator v1.7.4 [47]. Ten percent of the trees were discarded as burn-in, and the remaining trees were summarized as a maximum clade credibility tree including average divergence times and their associated 95% high posterior densities (HPD). Trees were visualized using FigTree v1.4.0 [46]. We used $K=r/2t$ to estimate the rate of substitutions per site [50].

## Demographic history

We used a Bayesian skyline plot [51] that was estimated by BEAST v1.7.4 [47] and mismatch distribution [52] to infer demographic history. Five independent runs with 30 million generations were conducted. The substitution model HKY+I with empirical base frequencies, was used with an uncorrelated lognormal relaxed clock model and a piecewise-constant coalescent Bayesian skyline tree prior with 10 starting groups. Trees and
parameters were sampled every 1000 iterations, with a burn-in of 10%. The results of each run were combined in LogCombiner v1.7.4 [47] and the result was visualized by using TRACER 1.5 [53]. In addition, mismatch distributions were obtained with the ARLEQUIN software package [38]. Mismatch distributions were calculated using the sudden expansion model [54] with 1000 parametric bootstraps. The sum of squared deviations (SSD) and Harpending's raggedness index (Hri) were calculated to assess the validity of the sudden expansion assumption.

## Results

### Sequence analysis, genetic diversity and differentiation

From the DNA samples of *M. gallopavo* collected for this study, we obtained 187 sequences (556 to 672 bp long) of the mtDNA D-loop that were registered in GenBank (Accession numbers: MF161996 to MF162182). Fifteen haplotypes were identified within the domesticated and wild turkey individuals analyzed in this study with an overall moderate Hd and low $\pi$ (Table 3), among which 11 were found in the domesticated turkeys from Mexico, Guatemala, and the commercial line Bronze with nine polymorphic sites, moderate Hd and low $\pi$. For the *M. g. mexicana* individuals, we detected five haplotypes with eight polymorphic sites, high Hd and low $\pi$. Finally, only one haplotype was identified in *M. g. intermedia*.

For the domesticated population from Mexico, we detected two dominant haplotypes designated MICH1 (n=107) and MICH2 (n=44) present in 61.49 and

25.28% of the individuals in the population, respectively. Interestingly, 95 domesticated individuals from Michoacan, eight from Puebla, all those from Guatemala, one from the commercial line Bronze, all *M. g. intermedia* individuals, and one *M. g. mexicana* individual, which was originally designated as carrying the MGMDgoB haplotype, shared the MICH1 haplotype. Thus, considering that MGMDgoB and MICH1 are the same haplotype or maternal line, it was hereafter designated as the MGMDgoB/MICH1 haplotype. In addition, 36 other domesticated individuals from Michoacan, two from Puebla, and six from the commercial line Bronze shared the MICH2 haplotype. These results revealed that many of the domesticated turkeys of Mexico and Guatemala and the individuals of the Bronze commercial line shared the same haplotypes; therefore, in the following analyses, they were treated as a single group called domesticated/commercial.

Next, to corroborate and strengthen our results, the mtDNA D-loop sequences of all domesticated (described as *M. gallopavo*) and wild individuals reported in the NCBI GenBank database were included in the following analyses (Table 2). The total population analyzed (n=612) showed overall high Hd and $\pi$ (Table 4). All the sequences of domesticated turkeys that were present in the NCBI GenBank database were included in the domesticated/commercial group. The analysis of domesticated/commercial turkeys showed moderate Hd and low $\pi$ (Table 4). Among the wild populations, diversity

<table><caption>**Table 3** Genetic diversity indices for the domesticated/commercial, *M. g. mexicana*, and *M. g. intermedia* turkeys included in this study</caption><thead><tr><th>Population</th><th>n</th><th>nt</th><th>H</th><th>Hd</th><th>π</th><th>S</th></tr></thead><tbody><tr><td>Domesticated/commercial</td><td>174</td><td>549</td><td>11</td><td>0.558</td><td>0.00153</td><td>9</td></tr><tr><td><i>M. g. mexicana</i></td><td>9</td><td>632</td><td>5</td><td>0.806</td><td>0.00483</td><td>8</td></tr><tr><td><i>M. g. intermedia</i></td><td>4</td><td>637</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Total</td><td>187</td><td>535</td><td>15</td><td>0.586</td><td>0.00198</td><td>16</td></tr></tbody></table>

n number of individuals, nt number of characters considered in the matrix, H number of haplotypes, Hd haplotype diversity, π nucleotide diversity, S number of polymorphic sites

<table><caption>**Table 4** Genetic diversity indices for each population</caption><thead><tr><th>Population</th><th>n</th><th>H</th><th>S</th><th>Hd</th><th>π</th><th>D-Tajima</th></tr></thead><tbody><tr><td>Domesticated/commercial</td><td>194</td><td>13</td><td>12</td><td>0.569</td><td>0.00162</td><td>-1.60009 (NS, 0.10 &gt; P &gt; 0.05)</td></tr><tr><td><i>M. g. intermedia</i></td><td>66</td><td>13</td><td>13</td><td>0.851</td><td>0.00530</td><td>-0.43092 (NS, P &gt; 0.10)</td></tr><tr><td><i>M. g. merriami</i></td><td>217</td><td>18</td><td>15</td><td>0.628</td><td>0.00826</td><td>1.09874 (NS, P &gt; 0.10)</td></tr><tr><td><i>M. g. mexicana</i></td><td>36</td><td>6</td><td>9</td><td>0.308</td><td>0.00209</td><td>-1.72617 (NS, P &lt; 0.10)</td></tr><tr><td><i>M. g. osceola</i></td><td>29</td><td>8</td><td>8</td><td>0.813</td><td>0.00402</td><td>-0.41858 (NS, P &gt; 0.10)</td></tr><tr><td><i>M. g. silvestris</i></td><td>70</td><td>22</td><td>23</td><td>0.913</td><td>0.00485</td><td>-1.71251 (NS, 0.10 &gt; P &gt; 0.05)</td></tr><tr><td>All samples</td><td>612</td><td>64</td><td>46</td><td>0.888</td><td>0.00774</td><td>-1.34597 (NS, P &gt; 0.10)</td></tr></tbody></table>

n number of individuals, H number of haplotypes, S polymorphic sites, Hd haplotype diversity, π nucleotide diversity
levels varied with the *M. g mexicana* population showing the lowest Hd and π levels, and the *M. g. silvestris* and *M. g. merriami* populations the highest levels (Table 4).

Analysis of the genetic differentiation between the domesticated/commercial group and the wild populations showed that there is differentiation among these populations (Table 5). The highest differentiation was observed between domesticated/commercial turkeys and *M. g. mexicana*. In addition, the lowest genetic differentiation values were found between *M. g. intermedia*, *M. g. osceola* and *M. g. silvestris*. Based on these results, we defined three groups of *M. gallopavo*: mexicana, merriami, and intermedia/silvestris/osceola/domesticated/commercial. The distribution of genetic variation obtained by AMOVA without a priori defined groups revealed that the genetic variation was highest within populations (Table 6). The percentage of genetic variation among populations and the high fixation index indicated a structure with subpopulations within species. When we divided the population into three groups, the percentage of genetic variation between them was equal to 31.14%. The fixation index $F_{CT}$ reached a high value indicating a high level of genetic differentiation among groups (Table 6).

## Haplotype network

We constructed a haplotype network to visualize the relationships between haplotypes and their frequencies for the 612 domesticated/commercial and wild turkey sequences (see Additional file 1: Table S1). The analyses revealed 64 haplotypes that differed from each other by a small number of mutations. The network (Fig. 2) shows eight haplogroups, each with a dominant haplotype. Four mexicana haplotypes clustered with those obtained from the NCBI GenBank database, forming a haplogroup with 34 individuals that shared the dominant haplotype Mgm (Fig. 2). A haplogroup that derives from haplotype Mgm has the dominant haplotype Mg, which is shared by intermedia, merrami, and silvestris samples. Merriami turkeys integrate a haplogroup with the dominant haplotype Mgmer that is shared with the archeological samples [7]. Haplotypes of the intermedia individuals are dispersed and shared with domesticated/commercial, osceola, and silvestris turkeys. A haplogroup that contains mainly silvestris turkeys showed a dominant haplotype (Mgs) that is shared with domesticated (from Puebla, Mexico), osceola, intermedia, and merriami individuals, and with peripheral haplotypes forming a star. Moreover, Mgs is related to the dominant haplotypes Mgo, MGM-DgoB/MICH1, and MgArch through one mutation. Two haplogroups were linked with the dominant haplotypes

<table><caption>**Table 5 Pairwise genetic differentiation ($F_{ST}$) of populations**</caption><thead><tr><th></th><th>M. g. merriami</th><th>Domesticated/commercial</th><th>M. g. silvestris</th><th>M. g. intermedia</th><th>M. g. mexicana</th><th>M. g. osceola</th></tr></thead><tbody><tr><td>M. g. merriami</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>Domesticated/commercial</td><td>0.50550</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>M. g. silvestris</td><td>0.25394</td><td>0.49995</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>M. g. intermedia</td><td>0.23189</td><td>0.35257</td><td>0.14670</td><td>–</td><td>–</td><td>–</td></tr><tr><td>M. g. mexicana</td><td>0.43530</td><td>0.85891</td><td>0.65913</td><td>0.56566</td><td>–</td><td>–</td></tr><tr><td>M. g. osceola</td><td>0.25278</td><td>0.54871</td><td>0.12507</td><td>0.09326</td><td>0.68575</td><td>–</td></tr></tbody></table>

<table><caption>**Table 6 Summary AMOVA without a priori defined groups and in three groups**</caption><thead><tr><th>Source of variation</th><th>d.f.</th><th>Sum of squares</th><th>Variance components</th><th>Percentage of variation</th><th>Fixation index</th></tr></thead><tbody><tr><td colspan="6"><em>Without a priori defined groups</em></td></tr><tr><td>Among population</td><td>5</td><td>387.469</td><td>0.83993</td><td>43.98</td><td></td></tr><tr><td>Within population</td><td>606</td><td>648.225</td><td>1.06968</td><td>56.02</td><td>F<sub>ST</sub>=0.43984</td></tr><tr><td>Total</td><td>611</td><td>1035.694</td><td>1.90961</td><td>100</td><td></td></tr><tr><td colspan="6"><em>In three groups<sup>a</sup></em></td></tr><tr><td>Among groups</td><td>2</td><td>299.146</td><td>0.65399</td><td>31.14</td><td>F<sub>CT</sub>=0.31139</td></tr><tr><td>Among populations within groups</td><td>3</td><td>88.323</td><td>0.37655</td><td>17.93</td><td>F<sub>SC</sub>=0.26036</td></tr><tr><td>Within populations</td><td>606</td><td>648.225</td><td>1.06968</td><td>50.93</td><td>F<sub>ST</sub>=0.49068</td></tr><tr><td>Total</td><td>611</td><td>1035.694</td><td>2.10021</td><td>100</td><td></td></tr></tbody></table>

<footnote><sup>a</sup> In three groups: (1) mexicana, (2), merriami, and (3) intermedia/osceola/silvestris/domesticated/commercial</footnote>
<figure><img src="image_4.png" /></figure>

MGMDgoB/MICH1 and MICH2. Haplotype MGMDgoB/MICH1 is shared by domesticated turkeys from Mexico, Guatemala, Canada, and the USA, individuals of the commercial line Bronze and wild turkeys in the intermedia and mexicana populations. It was also detected in individuals that were previously identified as wild *M. g. gallopavo*, for which samples were collected in 1903 in Veracruz and Michoacan, Mexico [7]. Additional haplotypes derived from MGMDgoB/MICH1 are present in domesticated turkeys of Michoacan, individuals of the commercial line Bronze, and mexicana, intermedia, and silvestris turkeys. Haplotype MICH2 is shared between domesticated turkeys from Mexico and Canada, most individuals of the commercial line Bronze, an osceola individual, and it corresponds to the same haplotype of an individual that was identified as wild *M. g. gallopavo* from Michoacan and collected in 1903 [7]. Derived from the MICH2 haplotype, peripheral haplotypes were identified that are present in the commercial line Bronze and domesticated turkeys from Michoacan (Fig. 2).

## Genealogical relationships

Phylogenetic analyses were performed to estimate the genealogical relationships between these haplogroups using the haplotypes that were detected in this study and the sequences that are designated as subspecies or domesticated in the NCBI GenBank database, which amounts to 64 haplotypes (Table 2). In the ML and BI consensus tree, haplotypes from the wild mexicana population, which was identified as a haplogroup in the network, are basal (C1) (Figs. 2, 3). Although the topology of the consensus tree shows polytomies, we can observe a clade (C2) that integrates intermedia and merriami haplotypes, which were detected as a haplogroup in the network, and a subclade that contains merriami haplotypes, which are linked with the haplogroup with the dominant haplotype Mgmer (Figs. 2, 3). In addition, a large clade (C3) includes haplotypes that were identified in domesticated/commercial, archeological, intermedia, merriami, silvestris, oceola, and mexicana (MGMDgoB/MICH1 and MGMDgoD) individuals, which belong
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Genealogy of *M. gallopavo* obtained with Bayesian inference (BI) and maximum likelihood (ML) analyses. Estimates were based on 441 characters of *M. gallopavo* mitochondrial control region DNA sequences. The out-group is the sister species *M. ocellata*. Names of distinct clades/groups are indicated. Values over the branches represent posterior probabilities and bootstrap values (PP/BP). (*) Values below PP = 0.5 or PB = 50. The keys of the taxa are in Table 2</figcaption></figure>

to five haplogroups detected in the haplotype network (Figs. 2, 3). In one polytomy, we found haplotypes of the Mgo, Mgs and MGMDgoB/MICH1 haplogroups. The C3 clade also contains six subclades or expansions: three major ones (SCI, SCII, and SCIII) and three with

two haplotypes each (SCIV, SCV, and SCVI). In SCI, the haplotypes corresponding to the haplogroup with the dominant haplotype MICH2 were included (MICH2, MICH5aqui, MICH8coah, Com9, MG80 and MG52) (Figs. 2, 3). SCII, which covers mainly archaeological
samples and one merriami individual, corresponds to the
haplogroup with the dominant haplotype MgArch. SCIII
is comprised of silvestris haplotypes, which correspond
to the expansion of the haplogroup with the dominant
haplotype Mgs (Figs. 2, 3). Finally, the three subclades
with two haplotypes each also correspond to expansions
of haplotypes shared by silvestris and osceola turkeys,
which in the network are related to the haplotype Mgs
through one mutation.

Analysis of divergence times

Using Bayesian inference, we analyzed divergence times
to estimate when the separation of groups and expan-
sions occurred. The results placed the most recent com-
mon ancestor (MRCA) of the genera *Meleagris* and
*Gallus* at 33.66 million years ago (Mya). *M. gallopavo* and
*M. ocellata* shared the MRCA during Pliocene-Miocene
time limits [5.35 Mya, HPD (95%) = 2.72–10.17]. Within
*M. gallopavo*, the differentiation of the basal group from
northwestern Mexico (mexicana population) began dur-
ing the Pliocene (3.39 Mya) (Fig. 3). The more ancestral
haplotypes of the intermedia and merriami turkeys share
a MRCA in the Pleistocene (1.65 Mya) (Fig. 3). In clade
C2, the subclade, which is composed exclusively of mer-
riami turkeys, originated 1.02 Mya during the Pleistocene
and underwent a subsequent diversification with two
haplotypes 0.43 Mya during the Pleistocene (Fig. 3). In
clade C3, diversification occurred at different times dur-
ing the Pleistocene. Subclades SCI with haplotypes that
were detected in domesticated turkeys originated 0.65
Mya, SCII with merriami and archeological haplotypes
diverged 0.39 Mya, and SCIII with silvestris haplotypes
originated 0.92 Mya (Fig. 3). The three other subclades
diverged 0.08 (SCIV), 0.28 (SCV), and 0.07 (SCVI) Mya

(Fig. 3). However, these results should be considered with
caution because of the low level of genetic variability in
the data analyzed.

Demographic history

Finally, to describe the changes in effective population
size through time, we investigated the demographic his-
tory of the species. The mismatch distribution graph
shows that the highest frequency of pairwise differences
is around 2, which indicates that the population analyzed
in this study has few mutations, and represents closely-
related individuals. The values recovered from SSD
(0.010, *P*=0.35) and Hri (0.027, *P*=0.37) indicate that
the analyzed data fit the sudden expansion model [54]
(Fig. 4a).

In the skyline plot, the y-axis shows population size expressed in N_e, where N_e is the effective population size and is the generation time, which is approximately one year in *M. gallopavo* [55, 56]. The time axis was scaled using the rate of 0.0046 substitutions/site/million years (SSM) that was obtained in this study. The skyline plot analysis of all *M. gallopavo* individuals as a population shows a demographic reconstruction starting from 300,000 years ago (Kya) (Fig. 4b). Between 300 to 130 Kya, the population remained stable, followed by period of growth and then again a stable period from 90 to 10 Kya. In addition, a slight population decline at approximately 10 Kya was observed in the skyline plot.

Discussion

Genetic diversity and differentiation

Based on the analysis of sequences generated for this
study and from the NCBI GenBank database, our results
show a high level of haplotype and nucleotide diversity

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Mismatch distribution and Bayesian skyline plots. Mismatch distribution for the domesticated group (a). Histograms correspond to the observed frequencies; the line represents expected frequencies under the sudden expansion model. Bayesian skyline plots show the population history of *M. gallopavo* (b), with the black line indicating median population size estimates expressed in N<sub>e</sub> through time; colored areas represent 95% HPD intervals</figcaption></figure>

<figure><img src="image_7.png" /></figure>
which suggests that the turkey population has remained stable with an old evolutionary history [57]. Analysis of metrics of genetic diversity by group revealed particu- lar histories. Domesticated/commercial turkeys showed moderate Hd and low π, which indicates that they origi- nated from a small number of founders [57]. Silvestris, osceola, and intermedia populations showed the highest Hd and a low π, which suggests a bottleneck followed by a rapid expansion [57], whereas the merriami popula- tion had a moderate Hd and a high π, which suggests that it has remained stable [57]. Regarding the sequences of mexicana individuals that were obtained from the NCBI GenBank database (n=27), all except one individual dis- played the Mgm haplotype that we detected in this work, which results in this population having both low Hd and π, suggesting a bottleneck (Table 4) [57]. Nevertheless, since we identified five haplotypes for nine mexicana individuals, we believe that increasing sample size and extending the geographic area analyzed would lead to the detection of additional haplotypes.

In agreement with the levels of pairwise genetic dif-
ferentiation and distribution of genetic variation by
AMOVA analysis, approximately three groups of *M. gal-
lopavo* were identified: (1) mexicana, (2) merriami, and
(3) intermedia/silvestris/osceola/domesticated/com-
mercial. This result indicates that although these groups
share haplotypes, the proportion of unshared haplotypes
exclusive to each group is significant and variable [58].

Haplotype network

We found 64 haplotypes that differ from each other by
a small number of mutations, with some of these haplo-
types being shared by domesticated and wild populations,
which indicates that in *M. gallopavo* there is no sub-spe-
ciation; this is in agreement with previous reports for
the species [5, 7, 14, 15]. Eight haplogroups, each with a
dominant haplotype, were identified. These haplogroups
corresponded with the three groups that were identified
by the analyses of genetic differentiation and distribu-
tion of genetic variation as follows: (1) mexicana (haplo-
group with the dominant haplotype Mgm), (2) merriami
(haplogroup with the dominant haplotype Mgmer), and
(3) intermedia/silvestris/osceola/domesticated/commer-
cial (haplogroups with the dominant haplotypes Mgs,
Mgo, Mg, MGMDgoB/MICH1, MICH2, and MgArch)
(Fig. 2). Considering that the haplotypes identified in
domesticated/commercial turkeys are shared with those
in wild turkeys, we designated the third group as mexi-
cana/intermedia/silvestris/osceola. Figure 5 shows the
geographic distribution of the haplotypes, which agrees
with the three detected groups. It should be noted that
the detection of the dominant haplotype Mgs in domes-
ticated turkeys from Mexico and wild populations

indicates that its distribution ranges from the central to
northeastern and southeastern USA, and from northeast-
ern to central Mexico (Fig. 5). The fact that the dominant
haplotype Mgs and its peripheral haplotypes form a star,
is typical of population expansions from a small number
of founders [57]. This coincides with the genetic diver-
sity analysis of the silvestris population, for which high
Hd and π levels and a negative, although not significant,
Tajima's D value (Table 4) were found, which also suggest
expansion of the population [57].

On the one hand, our results indicate that the domi-
nant haplotypes MGMDgoB/MICH1 and MICH2 are the
main founding maternal lines of domesticated turkeys.
On the other hand, since these two maternal lines were
detected in individuals of the commercial line Bronze
and in domesticated individuals from Canada and the
USA, they probably constituted the basis of the current
highly selected commercial lines. Our results agree with
those of a study that evaluated the genetic diversity of dif-
ferent turkey populations by using single nucleotide poly-
morphisms (SNPs), i.e. that the heritage turkey varieties
Royal Palm and Narragansett and commercial popula-
tions derive from the wild mexicana population, that the
commercial lines share the same origin, and that possibly
specific haplotypes (nuclear DNA) were selected in the
modern domesticated turkey [14]. In addition, our find-
ing that some individuals sampled from Guatemala carry
the MGMDgoB/MICH1 haplotype indicates that there
has been an exchange between northern Mesoamerica
and the Maya cultural region, as proposed by Thornton
et al. [59]. To confirm this, it would be necessary to ana-
lyze individuals from Central America.

Genealogical relationships and divergence time

Our estimate of the date of divergence between the com-
mon ancestor of *Meleagris* and *Gallus* genera i.e. 33.66
Mya agrees with the estimates reported by Claramunt
and Cracraft [60]. The separation between *M. gallopavo*
and *M. ocellata* 5.35 Mya coincides with the Miocene–
Pliocene boundary, which indicates a deep divergence
between lineages. We identified different diversifica-
tion events in turkey particularly during the Pleistocene,
which coincides with reports on fossils dated between
0.3 and 2.6 Mya [49]. At the base of the tree, mexicana
turkeys (haplogroup Mgm), which currently inhabit
northwestern Mexico, are a relict population of *M. gal-
lopavo*. Similar to the results reported by Mock et al. [5]
and Speller et al. [7], we found that the mexicana group,
represented by the MGM60, MGMDgoA, MGMDgoC,
and MGMDgoE haplotypes, is ancestral (C1 in Fig. 3).
Our results suggest a diversification process in the mexi-
cana population during the Pliocene (3.39 Mya). We pro-
pose that the range of this population expanded towards
<figure><img src="image_8.png" /></figure>

the north in Arizona and New Mexico and towards the
center of the USA. A second group (C2) that comprises
the intermedia and merriami haplotypes (haplogroup
Mg) originated in the same area (Figs. 2, 3, 5). Currently,
merriami and mexicana populations occupy ponderosa
pine and pine-oak woodlands of the southwestern USA
and northern Mexico, respectively. Apparently, the great
deserts of North America, are an efficient geographic
barrier for the groups detected in the current analysis. In
western North America (southwest of USA and north-
west of Mexico), we identified some genetic disconti-
nuities that are associated with the Sonora-Mojave and
Chihuahua deserts, which we suppose have isolated and

promoted the divergence between mexicana and merri-
ami populations.

Based on our results on the establishment of several
haplogroups in the geographical space, genealogical rela-
tions and genetic differentiation (Figs. 3, 5 and Table 5),
we propose that the mexicana/intermedia/silvestris/
osceola genetic group (C3) expanded its range from the
center of the USA east to the Atlantic coast and to the
south through the Sierra Madre Oriental (SMOr) until it
reached the center of Mexico (Figs. 3, 5). Genetic discon-
tinuities have been identified within different vertebrate
and plant species in the eastern of USA [61]. Based on
our results, in *M. gallopavo* there is no genetic discontinuity in the region (Fig. 5).

However, the occurrence of shared haplotypes in individuals from locations in northeastern Mexico (intermedia) to Michoacan led us to propose that the Trans-Mexican Volcanic Belt (TMVB) may be responsible for the expansion of the range of the C3 group (Figs. 1, 5). We assume that the Mexican Plateau acted as a geographic barrier that limited the contact between populations of mexicana and the C3 group (northeastern-northwestern of Mexico) (Figs. 1, 5).

Finally, we found that the haplotypes present in domesticated turkeys originate from the genetic group C3, which includes the MGMDgoB/MICH1 haplotype and its peripheral haplotypes, and by expansion from the SCI, which includes the MICH2 haplotype and its peripheral haplotypes. Thus, domesticated turkeys do not originate from an extinct subspecies *M. g. gallopavo*.

The presence of the MGMDgoB/MICH1 and MGMDgoD haplotypes in the mexicana relict population, indicate that the mexicana population and the C3 group have been in contact, probably in central-northwestern of Mexico (Figs. 1, 2, 3). However, we did not identify mexicana relict haplotypes in domesticated turkeys, thus this contact between the mexicana population and the C3 group probably took place through the expansion of the species from the center to the northwest of Mexico. It is important to note that a previous analysis based on microsatellite markers showed that the MGM-DgoB/MICH1, and MGMDgoD haplotypes are present in the wild mexicana population [11]. In addition, the fact that the samples of wild turkeys collected in 1903 in the Michoacan and Veracruz areas of Mexico shared the haplotypes MGMDgoB/MICH1 and MICH2 indicates that these old and new maternal lines persisted in the wild population of central Mexico until the last century. In contrast to other domesticated species for which events of interspecific hybridization or multiple origins have been observed [19, 26–28, 62, 63], the domestication of the turkey is less complex; according to our results, domestication of this species has a unique origin that is likely in the center of Mexico, whereas Nelson [10] and Leopold [6] both proposed that it was in Michoacan, Mexico. However, the low nucleotide diversity in the D-loop sequence of *M. gallopavo* makes it difficult to determine the center of origin with more precision.

## Demographic history

The rate of substitutions per site obtained in this study coincides with the range of substitutions per site per million years of mitochondrial genes for birds reported in the literature [64]. The multimodal-shaped mismatch

distribution suggests demographic stability (Fig. 4a). However, our data also support an expansion (SSD and Hri) that coincides with the population expansion observed in the skyline plot approximately 110 Kya (Fig. 4), which also coincides with the Eemian interglacial period (during the marine isotope stages MIS5e and MIS6d that occurred 133 to 103 Kya [65]), (Fig. 4). Our results of the analysis of genetic diversity (high nucleotide and haplotype diversities) support the observation that the population remained large and stable 90 to 10 Kya (Table 3). The observed population expansion followed by a stable period from 90 to 10 Kya is possibly associated with the expansions detected in the C3 group 80 and 70 Kya (Fig. 3, SCIV and SCVI).

The slight decline of the population about 10 Kya coincides with the cooling during the Younger Dryas. Moreover, it has been reported that 12,900 years ago an extraterrestrial impact occurred in northern North America that contributed to the late Pleistocene mega-faunal extinctions and adaptive shifts among Paleoamericans in North America [66]. Environmental changes caused by the combination of these events in North America may have impacted the availability of resources and consequently promoted the declines of the population observed in this study. Although it also is possible that human activities had some impact on the observed demographic decline, there is no evidence of intensive consumption of turkeys by Amerindian tribes (in USA), since they showed preference for large prey during hunting [67–69].

## Conclusions

A large and stable population of *M. gallopavo* occupied a wide geographical distribution from the north to the center of America (USA and Mexico) during the Pleistocene. Due to the expansion of their geographical range and to divergence events during the Pliocene–Pleistocene and Pleistocene–Holocene boundaries, three genetic groups originated within the species: mexicana, merriami and mexicana/intermedia/silvestris/osceola. The MGM-DgoB/MICH1 and MICH2 haplotypes and their peripheral haplotypes that belong to the mexicana/intermedia/silvestris/osceola group, are the main maternal lines that were captured for domestication in the center of Mexico (Trans-Mexican Volcanic Belt), the only region of turkey domestication. Domesticated turkeys populations from backyards in Michoacan come from the founding domesticated population. To confirm these results further, sampling of turkeys should be extended to key regions of Mexico (Sonora, Zacatecas, Jalisco, Nayarit, Tamaulipas and the center-south region). Finally, we provide new data on the haplotype diversity that prevails among
domesticated turkeys from backyards in Mexican rural
communities, with six haplotypes, which, to date, had not
been reported in *M. gallopavo*.

**Additional file**

**Additional file 1: Table S1.** Haplotypes identified, polymorphic sites and absolute frequencies for each population sampled, including sequences obtained from databases.

**Authors' contributions**

GPJ, RLZ and MEC carried out the experiments. GPJ carried out the bioinformatics analyses. MGZP, GPJ and HCC conceived and designed the study and drafted the manuscript. All authors read and approved the final manuscript.

**Author details**

¹ Centro Multidisciplinario de Estudios en Biotecnología, FMVZ, Universidad Michoacana de San Nicolás de Hidalgo, Km. 9.5 Carretera Morelia-Zinapécuaro, Posta Veterinaria, C.P. 58000 Tarímbaro, Michoacán, Mexico. ² Facultad de Medicina Veterinaria y Zootecnia, Universidad Autónoma de Tamaulipas, Km. 5, Carretera a Cd. Mante S/S, C.P. 87276 Ciudad Victoria, Tamaulipas, Mexico.

**Acknowledgements**

The authors thank the financial support provided by the Coordinación de la Investigación Científica, Universidad Michoacana de San Nicolás de Hidalgo (project to MGZP), and CONACYT (Projects CONACYT-SAGARPA 2004-C01-201, Fondos Mixtos CONACYT Gobierno del Estado de Michoacán 2009-05-115938, and scholarships granted to GPJ [No. 359650] and MECP [No. 226623]). The authors thank Dr. Omar Chassin Noria for his contribution with samples of domesticated turkeys from Guatemala.

**Competing interests**

The authors declare that they have no competing interests.

**Ethics approval and consent to participate**

Not applicable.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 31 May 2017 Accepted: 2 April 2018

Published online: 17 April 2018

**References**

1. Howell SN, Webb S. A guide to the birds of Mexico and northern Central America. Oxford: Oxford University Press; 1995. p. 225-6.

2. Porter R, Kirwan G. Wild turkey (*Meleagris gallopavo*). In: del Hoyo J, Elliott A, Sargatal J, Christie DA, de Juana E, editors. Handbook of the birds of the world alive. Barcelona: Lynx Edicions; 2017. http://www.hbw.com/node/53318. Accessed 25 Oct 2017.

3. Schorger AW. The wild turkey: its history and domestication. Oklahoma: University of Oklahoma Press; 1966.

4. Stangel PW, Leberg PL, Smith JI. Systematics and population genetics. In: Dickson JG, editor. The wild turkey: biology and management. Pennsylvania: Stackpole Books; 1992. p. 18-28.

5. Mock KE, Theimer TC, Rhodes OE Jr, Greenberg DL, Keim P. Genetic variation across the historical range of the wild turkey (*Meleagris gallopavo*). Mol Ecol. 2002;11:643-57.

6. Leopold AS. Wildlife of Mexico: the game birds and mammals. Berkeley: University of California; 1972. p. 268-75.

7. Speller CF, Kemp BM, Wyatt SD, Monroe C, Lipe WD, Arndt UM, et al. Ancient mitochondrial DNA analysis reveals complexity of indigenous

North American turkey domestication. Proc Natl Acad Sci USA. 2010;107:2807-12.

8. Monteagudo LV, Avellanet R, Azon R, Tejedor MT. Mitochondrial DNA analysis in two heritage European breeds confirms Mesoamerican origin and low genetic variability of domestic turkey. Anim Genet. 2013;44:786.

9. Crawford R. Introduction to Europe and the diffusion of domesticated turkeys from the Americas. Arch Zootec. 1992;41:307-14.

10. Nelson E. A Winter expedition in southwestern Mexico. Natl Geogr Mag. 1904;15:341-56.

11. Lopez-Zavala R, Cano-Camacho H, Chassin-Noria O, Oyama K, Vazquez-Marrufo G, Zavala-Paramo MG. Genetic diversity and population structure of Mexican domesticated turkeys. Rev Mex Cienc Pecu. 2013;4:417-34.

12. Kamara D, Gyenai KG, Geng T, Hammad H, Smith EJ. Microsatellite marker-based genetic analysis of relatedness between commercial and heritage turkeys (*Meleagris gallopavo*). Poult Sci. 2007;86:46-9.

13. Aslam ML, Bastiaansen JW, Elferink MG, Megens HJ, Crooijmans RP, Blomberg LA, et al. Whole genome SNP discovery and analysis of genetic diversity in turkey (*Meleagris gallopavo*). BMC Genomics. 2012;13:391.

14. Szalanski AL, Church KE, Oates DW, Bischof R, Powers TO. Mitochondrial-DNA variation within and among wild turkey (*Meleagris gallopavo*) subspecies. Trans Nebraska Acad Sci. 2000;26:47-53.

15. Mock KE, Theimer TC, Wakeling BF, Rhodes JOE, Greenberg DL, Keim P. Verifying the origins of a reintroduced population of Gould's wild turkey. J Wildl Manag. 2001;65:871-9.

16. Liu YP, Wu GS, Yao YG, Miao YW, Luikart G, Baig M, et al. Multiple maternal origins of chickens: out of the Asian jungles. Mol Phylogenet Evol. 2006;38:12-9.

17. Fumihito A, Miyake T, Takada M, Shingu R, Endo T, Gojobori T, et al. Mono-phyletic origin and unique dispersal patterns of domestic fowls. Proc Natl Acad Sci USA. 1996;93:6792-5.

18. Kanginakudru S, Metta M, Jakati RD, Nagaraju J. Genetic evidence from Indian red jungle fowl corroborates multiple domestication of modern day chicken. BMC Evol Biol. 2008;8:174.

19. Miao YW, Peng MS, Wu GS, Ouyang YN, Yang ZY, Yu N, et al. Chicken domestication: an updated perspective based on mitochondrial genomes. Heredity. 2013;110:277-82.

20. Zhang TJ, Li HF, Chen KW, Chang H, Tang QP, Zhang JX. Genetic diversity and systematic evolution of Chinese domestic ducks along the Yangtze-Huai River. Biochem Genet. 2007;45:823-37.

21. Zhang Y, Yang C, Ting Z, Huang ZY, Chen CY, Li XY, et al. Analysis of the genetic diversity and origin of some chinese domestic duck breeds. J Integr Agric. 2014;13:849-57.

22. He DQ, Zhu Q, Chen SY, Wang HY, Liu YP, Yao YG. A homogenous nature of native Chinese duck matrilineal pool. BMC Evol Biol. 2008;8:298.

23. Li HF, Zhu WQ, Song WT, Shu JT, Han W, Chen KW. Origin and genetic diversity of Chinese domestic ducks. Mol Phylogenet Evol. 2010;57:634-40.

24. Decker JE, Pires JC, Conant GC, McKay SD, Heaton MP, Chen K, et al. Resolving the evolution of extant and extinct ruminants with high-throughput phylogenomics. Proc Natl Acad Sci USA. 2009;106:18644-9.

25. Gautier M, Laloë D, Moazami-Goudarzi K. Insights into the genetic history of French cattle from dense SNP data on 47 worldwide breeds. PLoS One. 2010;5:e13038.

26. Decker JE, McKay SD, Rolf MM, Kim J, Molina Alcalá A, Sonstegard TS, et al. Worldwide patterns of ancestry, divergence, and admixture in domesticated cattle. PLoS Genet. 2014;10:e1004254.

27. Larson G, Dobney K, Albarella U, Fang M, Matisoo-Smith E, Robins J, et al. Worldwide phylogeography of wild boar reveals multiple centers of pig domestication. Science. 2005;307:1618-21.

28. Giuffra E, Kijas JM, Amarger V, Carlborg Ö, Jeon JT, Andersson L. The origin of the domestic pig: independent domestication and subsequent introgression. Genetics. 2000;154:1785-91.

29. Avise JC. Phylogeography: the history and formation of species. Cambridge: Harvard University Press; 2000.

30. Dutton PH. Methods for collection and preservation of samples for sea turtle genetic studies: NOOA Technical Memorandum; NOAA's National Marine Fisheries Service; 1996. p. 17-24.

31. López-Zavala R, Cano-Camacho H, Monterrubio-Rico T, Chassin-Noria O, Aguilera-Reyes U, Zavala-Páramo MG. Morphological and productive
characteristics of guajolote (*Meleagris gallopavo*) raised in backyard
systems in Michoacán, México. Livest Res Rural Dev. 2008;20:68.

32. FitzSimmons NN. Male marine turtles: gene flow, philopatry and mating systems of the green turtle *Chelonia mydas*. Ph.D. Thesis, University of Queensland; 1997.

33. Sanger F, Nicklen S, Coulson AR. DNA sequencing with chain-terminating inhibitors. Proc Natl Acad Sci USA. 1977;74:5463-7.

34. Genes Codes-Software Sequencher. http://www.genecodes.com/ sequencher. Accessed 06 Mar 2015.

35. Müller J, Müller K, Neinhuis C, Quandt D. PhyDE-Phylogenetic data editor. http://www.phyde.de/index.html. Accessed 13 March 2015.

36. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451-2.

37. Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics. 1992;131:479-91.

38. Excoffier L, Laval G, Schneider S. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinform Online. 2005;1:47-50.

39. Bandelt HJ, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37-48.

40. Fluxus Technology Ltd. Network 5. http://www.fluxus-engineering.com/sharenet.htm. Accessed 25 Jan 2017.

41. Posada D. jModelTest: phylogenetic model averaging. Mol Biol Evol. 2008;25:1253-6.

42. Alfaro ME, Huelsenbeck JP. Comparative performance of Bayesian and AIC-based measures of phylogenetic model uncertainty. Syst Biol. 2006;55:89-96.

43. Hasegawa M, Kishino H, Yano TA. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985;22:160-74.

44. Stamatakis A. RAxML Version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312-3.

45. Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19:1572-4.

46. Rambaut A. FigTree v1.4.0. http://tree.bio.ed.ac.uk/software/figtree/. Accessed 16 Mar 2015.

47. Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.

48. Kingman JFC. The coalescent. Stoch Process Appl. 1982;13:235-48.

49. Fossilworks. Gateway to the Paleobiology Database. http://fossilworks.org/. Accessed 10 Apr 2016.

50. Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol. 1980;16:111-20.

51. Drummond AJ, Rambaut A, Shapiro B, Pybus OG. Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol. 2005;22:1185-92.

52. Rogers AR, Harpending H. Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992;9:552-69.

53. Rambaut A, Suchard M, Drummond A. Tracer, version 1.5. http://tree.bio.ed.ac.uk/software/tracer/. Accessed 25 May 2016.

54. Schneider S, Excoffier L. Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA. Genetics. 1999;152:1079-89.

55. Leopold AS. The nature of heritable wildness in turkeys. Condor. 1944;46:133-97.

56. Blankenship LH. Physiology. In: Dickson JG, editor. The wild turkey: biology and management. Pennsylvania: Stackpole Books; 1992. p. 84-6.

57. Grant WS, Bowen BW. Shallow population histories in deep evolutionary lineages of marine fishes: insights from sardines and anchovies and lessons for conservation. J Hered. 1998;89:415-26.

58. Bird CE, Karl SA, Mouse PE, Toonen RJ. Detecting and measuring genetic differentiation. In: Held C, Koenemann S, Schubart DC, editors. Phylogenography and population genetics in Crustacea. Florida: CRC Press; 2011. p. 31-55.

59. Thornton EK, Emery KF. The uncertain origins of Mesoamerican turkey domestication. J Archaeol Method Theory. 2015;24:328-51.

60. Claramunt S, Cracraft J. A new time tree reveals Earth history's imprint on the evolution of modern birds. Sci Adv. 2015;1:e1501005.

61. Soltis DE, Morris AB, McLachlan JS, Manos PS, Soltis PS. Comparative phylogeography of unglaciated eastern North America. Mol Ecol. 2006;15:4261-93.

62. Jansen T, Forster P, Levine MA, Oelke H, Hurles M, Renfrew C, et al. Mitochondrial DNA and the origins of the domestic horse. Proc Natl Acad Sci USA. 2002;99:10905-10.

63. Luikart G, Gielly L, Excoffier L, Vigne JD, Bouvet J, Taberlet P. Multiple maternal origins and weak phylogeographic structure in domestic goats. Proc Natl Acad Sci USA. 2001;98:5927-32.

64. Pacheco MA, Battistuzzi FU, Lentino M, Aguilar RF, Kumar S, Escalante AA. Evolution of modern birds revealed by mitogenomics: timing the radiation and origin of major orders. Mol Biol Evol. 2011;28:1927-42.

65. Shackleton NJ, Sánchez-Goñi MF, Pailler D, Lancelot Y. Marine isotope sub-stage 5e and the Eemian interglacial. Glob Planet Change. 2003;36:151-5.

66. Firestone RB, West A, Kennett JP, Becker L, Bunch TE, Revay ZS, et al. Evidence for an extraterrestrial impact 12,900 years ago that contributed to the megafaunal extinctions and the Younger Dryas cooling. Proc Natl Acad Sci USA. 2007;104:16016-21.

67. Hill ME. Variation in Paleozoian fauna use on the Great Plains and Rocky Mountains of North America. Quat Int. 2008;191:34-52.

68. Hill ME Jr. A moveable feast: variation in faunal resource use among central and western North American Paleozoian sites. Am Antiq. 2007;72:417-38.

69. Lyman RL. Paleozoian exploitation of mammals in eastern Washington State. Am Antiq. 2013;78:227-47.

70. Crandall KA, Templeton AR. Empirical tests of some predictions from coalescent theory with applications to intraspecific phylogeny reconstruction. Genetics. 1993;134:959-69.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_9.png" /></figure>