Phylogeographic and phenotypic divergence between two subspecies of Testudo graeca (T. g. buxtoni and T. g. zarudnyi) across their contact zone in Iran

Contact zones are considered as windows into the evolutionary process, allowing identification of factors influencing the evolutionary forces. Here, we combined phylogenetic and morphometric analyses to explore the evolutionary process affecting the taxonomic pattern of two subspecies of Testudo graeca (T. g. buxtoni and T. g. zarudnyi) across their contact zone in Central Iran. Our results showed high levels of phylogeographic and phenotypic variation in the contact zone. Two monophyletic clades including, clade 1 (T. g. zarudnyi) and clade 2 (T. g. buxtoni) were identified. Furthermore, four distinct subclades were found in T. g. buxtoni, across a wide geographic range. Divergence time analysis suggests that the two subspecies diverged from one another after the uplifting of the Zagros Mountains during the early Pliocene. Using neutrality tests and mismatch distribution analysis, we found no evidence of recent population expansion. Morphological associations among geographical populations in the contact zone found more distinctions, with some significant adaptive and non-adaptive morphological variations in these populations. These distinctive morphological populations can be considered as management units (MUs) to conserve the evolutionary potential of this species. Finer scale evolutionary studies are required to address the southern part of the Zagros mountain range, where the overlapping of mitochondrial clades and subclades has occurred. Such information is essential for effective conservation of T. graeca populations, preventing translocation or mixing of individuals without comprehensive genetic and morphological assessment.


Results
Mitochondrial DNA: genetic diversity and phylogeographic analyses. Phylogenetic status of Iranian spur-thighed tortoises from the contact zone between T. g. buxtoni and T. g. zarudnyi was evaluated, relative to other subspecies of the species' global range ( Fig. 1), by sequencing a 948 bp fragment of the mitochondrial cytochrome b gene. In total, 37 sequences of the two Iranian subspecies were obtained. These sequences were aligned with 229 published spur-thighed tortoise cytochrome b sequences, retrieved from GenBank. See Table S1 for a complete list of sequences and corresponding haplotypes used in the present study. All new sequences were submitted to GenBank (accession numbers: ON887288-ON887307).
Based on 948 base pairs of cytochrome b from 336 mtDNA sequences, 128 haplotypes, including five new haplotypes from Iran, were identified. Analyses of maximum likelihood (ML) and Bayesian inference (BI) yielded identical tree topologies in which 10 lineages from the species' worldwide range (Fig. 2) were recognized. Spurthighed tortoises were divided into eastern Mediterranean (clades 1 to 5) and western Mediterranean (clades 6 to 10) lineages.
In the contact zone, two monophyletic clades: clade 1 (subspecies T. g. zarudnyi), and clade 2 (subspecies T. g. buxtoni) were identified. Clade 1 is predominantly distributed throughout the eastern part of the contact zone, whereas clade 2 is mainly distributed in the western part of the contact zone with an overlap across the central part of the contact zone (Fig. 1). Clade 2 is further subdivided into four subclades (A, B, C, and D, Fig. 2), three of which (A, C, and D) were found in Iran, whereas subclade B was only identified in sequences from Turkey. The subclade C is the most geographically wide-spread group in the western part of the contact zone, however, the other two subclades (A and D) are restricted to small geographic areas ( Fig. 1) (Table 1).
In an analysis of molecular variance (AMOVA), significant F ST was obtained among populations, with 92% of the variation was partitioned among populations and 8% within populations. The significant genetic structure was proved by high fixation index (F st = 0.92) ( Table 2).

Molecular clock estimates.
Divergence times of all supported nodes (> 0.5/50) in the phylogenetic tree are given in Table 3. The divergence time for T. graeca in Iran was averaged at about 6.7 Mya with 95% highest posterior density (HPD) intervals of 5-8.4. The divergence time of T. g. zarudnyi was estimated at about 3.9 Mya (95% HPD: 2.3-5.5), while the value for T. g. buxtoni was about 4 Mya (95% HPD: 2.7-6.9).
Demographic analysis and neutrality tests. Mismatch distribution analyses and neutrality tests suggested no evidence for unimodal distribution and didn't reveal significant population expansion (Fig. 4).
Traditional morphometric analyses. The cluster analysis revealed the presence of four distinct morphological clusters across the study area (Fig. 5). When the number of clusters was changed, the degree of similarity significantly decreased and the distance increased. Therefore, the optimal number of clusters was four ( Table 4). The sampling localities were grouped into four populations, including Central part of Zagros Eco-Region (CZ), Southern part of Zagros Eco-Region (SZ), Plain areas of Irano-Turanian Eco-Region (P.I-T) and Mountainous areas of Irano-Turanian Eco-Region (M.I-T) (Fig. 5).
Geometric morphometric analyses. In order to assess the shape changes of the carapace and plastron, their TPS files were imported separately into PAST and Morpho J software. Three CVA components were extracted. The first and second components accounted for 79.2% and 77.4% of variation in carapace and plas-  Graciá et al. 10,28 , Turkozan et al. 9 , and Javanbakht et al. 7 . (B) Geographical distribution of the sampling sites of the spur-thighed tortoise in the contact zone between two subspecies T. g. zarudnyi (Clade 1) and T. g. buxtoni (Clade 2) in Central Iran. Purple squares and red circles indicate clades 1 and 2, respectively, and the numbers next to the symbols indicate the haplotypes detected in each sampling sites, according to the numbering assigned to them (see Figs. 2, 3 for details). A, C, and D, denote three subclades found in clade 2 from the study area. Outlined symbols show the sampling sites of the species in Javanbakht et al. 7  www.nature.com/scientificreports/ tron, respectively. The P values (P < 0.0001) from permutation tests (10,000 permutation rounds) indicated significant Mahalanobis distances among studied populations (Table 6). In plastron shape data, overlap was observed between M.I-T and CZ populations, while no overlap was observed in CZ and P.I-T (Fig. 7a). Based on the carapace data, only M.I-T and SZ populations showed overlap (Fig. 7b). The most distinct population based      www.nature.com/scientificreports/ on carapace shape variation was CZ, while the maximum Mahalonobis distance was observed between P.I-T and other populations, based on the plastron shape variation. To further illustrate shape changes, transformation grids and wireframe graphs were applied on the carapace (Fig. 8) and plastron ( Fig. 9) data. In the first component, the main carapace shape changes were in the supracaudal and fifth vertebral scales, while NL was considerably changed in the second component. The middle landmarks showed the least amount of displacement. The highest displacement, in the first plastron CVA component, was related to landmark No. 18 (connection of femoral and anal sutures). In addition to this landmark, the most movements corresponded to landmark No. 4 in the second component. Cluster dendrograms were constructed based on the consensus of each studied population. The Copernic coefficient of this analysis for   www.nature.com/scientificreports/   www.nature.com/scientificreports/ carapace and plastron were 0.92 and 0.87 respectively, suggesting the distinction of studied populations based on both carapace (Fig. 10a) and plastron (Fig. 10b).

Discussion
The phylogeography, morphological diversity, and population structure of the most widely distributed tortoise, T. graeca have been investigated by several studies 3,4,6-9,18,29 . However, contact zones among subspecies can raise questions regarding the evolutionary boundaries and admixture 28 . Therefore, the current study concentrated  www.nature.com/scientificreports/ on the contact zone of two early-divergent and endemic subspecies of T. graeca (T. g. buxtoni and T. g. zarudnyi) in Central Iran to fill phylogeographic and phenotypic gaps in this area. Moreover, the results provide essential information to develop effective conservation plans. The main results are elaborated in the sections below.
Mitochondrial diversity and phylogeography. Our finding confirmed the general topology suggested by previous studies 3, 10 . Incorporating sufficient samples of eastern subspecies into our analyses has led to some achievements in the phylogenetic structure of the species. In our analysis, 10 divergent clades were identified, associating with the current known subspecies of T. graeca. The two subspecies, T. g. zarudnyi, and T. g. buxtoni formed two monophyletic clades. Four subclades were further identified within T. g. buxtoni, three of which (A, C, and D) were found in Iran, whereas subclade B was only identified in sequences from Turkey. The subclade C was the most geographically wide-spread group in the contact zone. However, the other two subclades (A and D) were restricted to small geographic areas. Subclades A and B could probably adapt to their local environment, while subclades C and D (also defined as haplogroups in previous studies 7 ) were widely distributed throughout T. g. buxtoni geographic range. Subclade D was mainly distributed outside the contact zone, with a high level of genetic diversity. Subclade C was the most diverse and widely distributed haplotype in the contact zone. It shares one haplotype (H94) in the northern and southern parts of T. g. zarudnyi distribution range. We found T. g. buxtoni, the only subspecies, with distinct subclades in a wide geographic range. Therefore, it is recommended to introduce new subspecies, particularly in western part of the region, which has been pronounced in a considerable genetic diversity and differentiation through a long evolutionary history 9 .
The second subspecies (T. g. zarudnyi) had a smaller distribution range in the contact zone, having three haplotypes. Haplotype H5 is distributed further from the two subspecies border and confined to the previously recognized distribution range 7 , while haplotype H3 is the most frequent haplotype in the contact zone. This haplotype was also recorded in some localities in the distribution range of T. g. buxtoni, close to the border of the two subspecies.
Divergence time analysis suggests that the two subspecies diverged from one another after the uplifting of the Zagros Mountains during the early Pliocene 30 . The divergence time estimates were coherent with the study of Fritz et al. 3,4 and Garcia 10 , suggesting that T. graeca had spread over the Pliocene, though our dating results were in narrower ranges compared to estimates of the previous studies 10 .
Our results did not support the hypothesis of population expansion of T. graeca in the study area, which is in line with the previous study 7 . Therefore, the high genetic diversity as well as divergence, especially in the contact zone, could be the result of long term stability in their distribution ranges 7 .
With our dataset in the contact zone, we obtained high level of genetic diversity and divergence with some expected introgressions. Genetic diversity statistics indicate high values of haplotype and nucleotide diversity, particularly for T. g. Buxtoni (Table 1). Clear genetic structure depicted in the phylogenetic tree (Fig. 2) and the haplotype network (Fig. 3), particularly for T. g. zarudnyi, was confirmed by F st and an analysis of variance (Table 3). Our analyses showed the presence of divergent clades with no evidence of recent population expansion and frequent and widespread haplotypes (e.g. H94 and H3 as potential ancestral haplotypes) in the contact zone. Therefore, we recognized the contact zone of T. g. buxtoni and T. g. zarudnyi as an evolutionary significant zone in the distribution range of the species.
Traditional and geometric morphometrics. Our morphological study is the first to be performed in the contact zone between the two subspecies of T. graeca in Iran. We examined a number of morphological variables and identified the most important morphological features pronounced in shapes. We identified four distinct groups in Central part of Zagros Eco-Region (CZ), Southern part of Zagros Eco-Region (SZ), Plain areas of Irano-Turanian Eco-Region (P.I-T) and Mountainous areas of Irano-Turanian Eco-Region (M.I-T). These regions are different based on the most significant climatic factors, including annual precipitation and temperature seasonality (Table S2).
The high discriminatory power of geometric morphometrics 31,32 was led to differentiate morphological populations, using carapace or plastron landmarks. Our finding is in agreement with the suggestions of population-specific variation, phenotypic and ecological plasticity among populations of tortoises 4,6,33-35 . The greatest morphological distance was found between the Zagros (mainly T. g. buxtoni) and Irano-Turanian (mainly T. g. zarudnyi) populations. In geometric morphometric analyses, the differentiation among morphological populations was completely clear for CZ in carapace and P.I-T in plastron. Although Copernic coefficients represented the distinctiveness of morphological populations, based on the existing overlaps, some similarities between CZ and M.I-T populations (in plastron) and between SZ and M.I-T (in carapace) were observed. Moreover, SZ was located at the center of the traditional CVA scatter plot (Fig. 6), and in the same group with T. g. zarudnyi populations in the cluster analysis (Fig. 10). This implies that some parts of the SZ may have facilitated gene flow among tortoise populations, due to the lack of biogeographical barriers such as mountain ranges and rivers. For example, our cluster analyses suggest CZ and SZ are considerably distinct, possibly due to the presence of the Khersan river. In addition, different environmental conditions in plain and mountainous areas seem to influence the distinction between P. I-T and M. I-T (especially in plastron shape changes). For example, our genetic assessments demonstrate no difference between two morphological populations, P. I-T and M. I-T. This can be explained by the role of phenotypic plasticity in adaptation of the populations to novel environments 4,36-38 .
As it can be inferred from our results, environmental conditions act effectively on costal scales; carapace height, and width of plastron lobes 17,39 . In open areas with limited shelter and food, as is the case in P. I-T, digging ability is positively selected 40 as a more flatted carapace may improve burrowing. Further, the carapace (CW and particularly MCW) and plastral lobes were narrower in P. I-T populations to provide wider openings www.nature.com/scientificreports/ for head and limbs and facilitate mobility. By contrast, the height of the carapace, width of vertebral scales and plastral lobes are greater in the populations of CZ and M. I-T and they have longer coastal scales. This can be attributed to the different environmental conditions in CZ and M. I-T compared to P. I-T (e.g. more precipitation and lower temperature seasonality), which can stimulate more plant growth, and the availability of more shelters in mountainous landscapes. Our results illustrated the importance of supracaudal and nuchal carapace scales for distinguishing populations. Previous studies 39,41 have also recommended the supracaudal and nuchal scales as reliable diagnostic characters for distinguishing the populations of Testudo species. Furthermore, some morphological variables such as MGSL, AbSL, FSL, GSL, PSL and IHASO can be used for discriminating the populations. Among these variables, IHASO is visually more obvious. The decrease in inner height of the anterior shell in P.I-T and M.I-T is evidently due to the inclination of the gular scales inward.
Combining mitochondrial and morphological data. The separation of the two subspecies is remarkably proved by mitochondrial clades as well as traditional and carapace geometric morphometrics. It seems that the Zagros Mountains act as ecological barriers, separating the geographic ranges of the two subspecies, though geographical overlap between mitochondrial clades and morphological populations in the southern part of Zagros implies the possibilities of gene flow in this part of the contact zone. For example, similar environmental conditions of M. I-T and CZ may have resulted in similar morphological features in the populations of these two subspecies (see above). Based on the current study, we recommend that these distinctive morphological populations can be considered as management units (MUs) to conserve the evolutionary potential causing adaptive and non-adaptive morphological variations, more effectively 42 . These key phylogeographic and phenotypic achievements remind the importance of contact zones as suitable sites for studying evolutionary processes. However, finer scale evolutionary studies are required to address the southern part of the Zagros mountain range, where the overlapping of mitochondrial clades and subclades as well as the morphological populations are happening. Consequently, translocation or mixing of individuals, even for conservational purposes, cannot be acceptable without comprehensive genetic and morphological assessment.

Materials and methods
Ethic statement. Samples and data were collected in accordance with the guidelines and regulations of the Iranian Department of Environment (DOE) under permission No. 97/6549. All individuals were sampled, photographed, morphologically measured and released at the point of capture with no harm. Sample collection followed the ethical recommendations in the ARRIVE guidelines 43 . Biopsy sampling technique as a nondestructive, quick and safe method was carried out to reduce the possibility of bleeding and infection 44 . However, alcohol was used before and after the sampling.
Sampling. The sampling was conducted in two consecutive years (2018 and 2019) during the species activity period in the contact zone of T. g. buxtoni and T. g. zarudnyi in Central Iran (Fig. 1). A total of 86 individuals were collected from Fars, Kerman, Isfahan, Yazd, Chaharmahal and Bakhtiari and Kohgiluyeh and Boyer-Ahmad provinces. For morphological analyses, however, only 82 adult individuals were measured and photographed. Growth rings on the costal scale of the carapace were counted and individual with more than 12 years old were considered as adults 45 . Tissue samples were taken from individuals covering the contact zone, using a biopsy punch. Small biopsy punches (2 mm) were used to obtain 0.5 mm deep tissues from the scales under feet close to the spurs. The tissue was immediately stored in pure 99% ethanol. Laboratory procedures. Genomic DNA was extracted from tissue samples using WizPrep™ gDNA Mini Kit (Cell/Tissue) following the manufacturers' instructions. Polymerase chain reaction was performed for amplification of a 948 bp fragment of the mitochondrial cyt b gene, using CytbG (5′ AAC CAT CGT TGT WAT CAA CTAC 3′) and mt-f-na (5′AGG GTG GAG TCT TCA GTT TTT GGT TTA CAA GAC CAA TG 3′) primers 46,47 . Amplifications were performed in 25 µl volumes, containing 10 µl of the Ampliqon PCR master mix, 1 µl of each primer, 9 µl H 2 O and 4 µl of DNA. The thermocycling was performed using an initial denaturation at 95° for 5 min followed by 35 cycles of 30 s at 94°, 50 s at 55°, and 1 min at 72° and a final extension at 72° for 7 min. Purified PCR products were sequenced on an ABI PRISM 3730xl using the BigDye Terminator Cycle Sequencing kit v.3.1 (Applied BioSystems).
Alignment and sequence analyses. Sequences were edited with SeqScape v.2.6 software (Applied Biosystems) and aligned with 229 previously published mtDNA sequences from GenBank using the ClustalW algorithm implemented in Mega v.7 48 and checked visually. In addition to T. g. buxtoni and T. g. zarudnyi sequences, all other subspecies sequences were retrieved from GenBank and included in the final dataset. Sequences of T. hermanni (AJ888362), T. marginata (AJ888333), and T. kleinmanni (AJ888371) were used as outgroup representatives.
The number of polymorphic sites (S), nucleotide diversity (π), number of haplotypes (H), haplotype diversity (h), and the average number of pairwise differences (K) were calculated using DnaSP 5.10.1 49 . The substitution saturation was checked in DAMBE v.6 50 .
PartitionFinder v.2.1.1 51 was used to reveal the optimal partitioning scheme and nucleotide substitution model, using the Bayesian Information Criterion (BIC). Nucleotide substitution models include K80 + G for the first, HKY + G for the second, and GTR + G for the third codon position. Maximum likelihood (ML) and Bayesian inference (BI) approaches were employed to assess the phylogenetic relationships. Bayesian phylogenetic tree was constructed in MrBayes v. 3 www.nature.com/scientificreports/ simultaneously for 10 million generations and sampling every 1,000 generations. The first 25% of the sampled trees and estimated parameters were discarded as burn-in.  54 for 10,000 bootstrap steps. Two methods for species delimitation were then used to identify the specific boundaries in the species. These approaches were implemented on online servers: the Poisson Tree Process model (PTP, https:// speci es.h-its. org/ ptp/) and the General Mixed Yule Coalescent method (GMYC, https:// speci es.h-its. org/ gmyc/).
A Median-Joining (MJ) network analysis was implemented in PopArt v.1.735 55 to clarify the phylogenetic relationships among haplotypes. Pairwise F st and an analysis of molecular variance (AMOVA) were performed in Arlequin v.3.5.2.2 56 with 10,000 permutations to estimate the level of genetic structure among clades and subclades in the contact zone.

Molecular clock analysis.
Divergence times of the lineages were estimated using a relaxed lognormal clock model in Beast v. 2.6.4 57 . Two calibration dates of 7.6 and 9 Mya were respectively applied for inferring the minimum ages of the lineages. The first one consisted of T. kleinmanni plus T. marginata and the second was located between T. hermanni and the first calibration node 10 (the two red circles on the phylogenetic tree, Fig. 2). The best fitting nucleotide substitution model was set for each codon position separately based on the optimal partitioning scheme as mentioned above. The speciation Yule Process (YP) was selected as a tree prior. The MCMC analyses were run for 10 million generations with sampling each 1000 generations and convergence diagnostics were assessed using Tracer v.1.7.1 53 .
Demographic analysis and neutrality test. Mismatch distributions were estimated separately for T. g. zarudnyi and T. g. buxtoni to test if their frequency graph shows a chaotic/multimodal pattern characteristic for populations in demographic equilibrium, or an unimodal profile which is found in populations that experienced recent geographic expansion. The test was performed under the null hypothesis that the observed data fit the sudden expansion model, using a generalized least-square approach with 1000 bootstrap replicates. Other statistics for analyzing population expansions or declines were also calculated, i.e., Fu's FS 58 , Tajima's D 59 and Ramos-Onsins and Rozas's R2 60 using DnaSP 5.10.1 49 . The negative value of these statistics is interpreted as indicating recent demographic expansion.

Traditional morphometrics. Body measurements (~ 40 measurements) on different scales of carapace
and plastron were performed, using calipers with an accuracy of 0.02 mm. Measured carapace scales consisted of nuchal, five vertebral, and supracaudal scales in the middle, and eight pleural scales on both sides of the vertebral. In plastrons, each scale is divided into two parts by a suture, and the scutes of the plastron include gular, humeral, pectoral, abdominal, femoral and anal.
Traditional morphometric measurements were performed following the protocol of Türkozan et al. 18  Since the morphological features may change with increasing body size, the effect of dimensions must be reduced and differences between the groups should be due to the differences in body shape and not the relative size. Thus, standardizing morphometric data is required to reduce the changes caused by allometric growth 61 . In order to eliminate the effect of size, morphometric data were standardized prior to performing analyses, using the following formula. www.nature.com/scientificreports/ A cluster analysis was performed to find the optimal number of morphological populations. The sampling localities were grouped into four areas (representative of different landscapes and climatic conditions), including Central part of Zagros Eco-Region (CZ), Southern part of Zagros Eco-Region (SZ), Plain areas of Irano-Turanian Eco-Region (P.I-T) and Mountainous areas of Irano-Turanian Eco-Region (M.I-T). We also recorded the most significant factors shaping the distribution pattern of the individuals, i.e., annual precipitation and temperature seasonality of each area 7,62 to assess the differences, using one-way analysis of variance (ANOVA) ( Table S2). Inter-group variance was evaluated using canonical variate analysis (CVA). Furthermore, a discriminant analysis was performed to classify individuals into known populations. All analyses were performed, using the software PAST 63 and MINITAB 18 64 . Geometric morphometrics. Tortoises were labeled and photographed alive without anesthesia, using a Canon 350 D digital camera (8 MP). Photographic conditions such as camera settings, magnification, focus, lens size, camera distance from the sample surface and backlight were consistent for all individuals and images were taken from ventral and dorsal views. The camera was placed parallel to the ground plane, as such that the focal axis of the camera was parallel to the horizontal plane of reference and centered on the plastron and carapace.
On each image, 2D coordinates of fixed landmarks were digitized by the same person, using TPS DIG v. 2.31 65 . In total, 31 and 25 two-dimensional (2D) landmarks were digitized on carapace and plastron, respectively (Fig. 11). Landmarks were digitized on the intersections of dermal scales. The distribution of landmarks covers nuchal, vertebral, pleural, and supracaudal scales of the carapace and twelve left and right scales of plastron (gular, humeral, pectoral, abdominal, femoral and anal). To minimize observer-induced errors, one person digitized the landmarks. Landmarks were, then, analyzed using generalized Procrustes analysis (GPA) 66,67 , which mathematically removes the effects of non-shape variation and normalizes shape data at equal scale, allowing for an accurate comparison of shapes regardless of their position, orientation and isometric size. To evaluate the differences among populations, CVA was performed on the shape data. Finally, cluster analysis (CA) was performed to group the studied populations. All analyses were conducted using PAST 63