Evolutionary history of Mexican domesticated and wild Meleagris gallopavo

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 Biotecnologia (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).

Country	Sample types	State	Locality	Sample size
Mexico	M. g. mexicana	Durango	Canatlan (UMA El Durangueño)	9
	M. g. intermedia	Tamaulipas	Villa de Casas	4
	Domesticated	Michoacan	Bajio Region (Penjamillo, Indaparapeo, Zacapu, La Piedad, Jimenez, Capula, Churintzio, Zamora)	33
			Trans-Mexican Volcanic Belt Region (Zitacuaro, Cuitzeo, Benito Juarez, Tuxpan, Tacambaro, Ario de Rosales, Salvador Escalante)	31
			Balsas Region (Churumuco, Huetamo, Tiquicheo, Tepalcatepec, Tuzantla, Nuevo Urecho, Buenavista, Apatzingan)	26
			Sierra Region (Tumbiscatio, Arteaga, Aguililla)	28
			Costa Region (Coahuayana, Aquila, Arteaga)	26
		Puebla	Rafael Lara Grajales (Benito Juarez Ejido)	1
			Totimehuacan (San Baltazar Azumiatla)	5
			Soltepec (San Mateo)	7
			Hueyotlipan (Sto. Tomas Hueyotlipan)	4
	Commercial line Bronze	Michoacan	Ario de Rosales	9
Guatemala	Domesticated	Izabal	Izabal	4
Total				187

DNA extraction and amplification

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 F_ST 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 PaleoDB fossil database (82,258) [49] were used. Markov chain Monte Carlo (MCMC) analyses were run for 10 million generations, 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.

Sequence analysis, genetic diversity and differentiation

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.

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, oceola, 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, MGMDgoB/MICH1, and MgArch through one mutation. Two haplogroups were linked with the dominant haplotypes 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 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 expansions occurred. The results placed the most recent common 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 during 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 merriami 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 during 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 history 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.

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 which suggests that the turkey population has remained stable with an old evolutionary history [57]. Analysis of metrics of genetic diversity by group revealed particular histories. Domesticated/commercial turkeys showed moderate Hd and low π, which indicates that they originated 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 population 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 displayed 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 differentiation and distribution of genetic variation by AMOVA analysis, approximately three groups of M. gallopavo were identified: (1) mexicana, (2) merriami, and (3) intermedia/silvestris/osceola/domesticated/commercial. 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 haplotypes being shared by domesticated and wild populations, which indicates that in M. gallopavo there is no sub-speciation; 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 distribution of genetic variation as follows: (1) mexicana (haplogroup with the dominant haplotype Mgm), (2) merriami (haplogroup with the dominant haplotype Mgmer), and (3) intermedia/silvestris/osceola/domesticated/commercial (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 mexicana/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 domesticated turkeys from Mexico and wild populations indicates that its distribution ranges from the central to northeastern and southeastern USA, and from northeastern 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 diversity 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 dominant 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 different turkey populations by using single nucleotide polymorphisms (SNPs), i.e. that the heritage turkey varieties Royal Palm and Narragansett and commercial populations 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 finding 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 analyze individuals from Central America.

Genealogical relationships and divergence time

Our estimate of the date of divergence between the common 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 diversification 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. gallopavo. 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 mexicana population during the Pliocene (3.39 Mya). We propose that the range of this population expanded towards 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 northwest of Mexico), we identified some genetic discontinuities that are associated with the Sonora-Mojave and Chihuahua deserts, which we suppose have isolated and promoted the divergence between mexicana and merriami populations.

Based on our results on the establishment of several haplogroups in the geographical space, genealogical relations 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 discontinuities 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 MGMDgoB/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 megafaunal 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].