Genetic and morphological divergence among three closely related Phrynocephalus species (Agamidae)

Background The Qinghai-Tibetan Plateau (QTP) is the world’s highest and largest plateau, but the role of its uplift in the evolution of species or biotas still remains poorly known. Toad-headed lizards of the reproductively bimodal genus Phrynocephalus are a clade of agamids, with all viviparous species restricted to the QTP and adjacent regions. The eastern part of the range of the viviparous taxa is occupied by three closely related but taxonomically controversial species, P. guinanensis, P. putjatia and P. vlangalii. Here, we combined genetic (mitochondrial ND4 gene and nine microsatellite loci), morphological (11 mensural and 11 meristic variables), and ecological (nine climatic variables) data to explore possible scenarios that may explain the discordance between genetic and morphological patterns, and to test whether morphological divergence is associated with local adaptation. Results We found weak genetic differentiation but pronounced morphological divergence, especially between P. guinanensis and P. vlangalii. Genetically, the species boundary was not so clear between any species pair. Morphologically, the species boundary was clear between P. guinanensis and P. vlangalii but not between other two species pairs. Body size and scale characters accounted best for morphological divergence between species. Morphological divergence was related to habitat types that differ climatically. Conclusions Our study provides evidence for genetic and morphological divergence among the three closely related viviparous species of Phrynocephalus lizards, and supports the idea that natural selection in spatially heterogeneous environments can lead to population divergence even in the presence of gene flow. Our study supports the hypothesis that the evolutionary divergence between viviparous Phrynocephalus species was a consequence of environmental change after the uplift of the QTP. Electronic supplementary material The online version of this article (10.1186/s12862-019-1443-y) contains supplementary material, which is available to authorized users.


Background
Genetic divergence and speciation can occur in different parts of an ancestral species' range and even within habitats [1]. Genetic divergence within and among species is not always accompanied by clear phenotypic (morphological, anatomical, physiological, and/ or behavioral) differences due to silent mutations or phenotypic convergence [2]. However, it can give rise to significant phenotypic changes due to novel adaptations via selection that drives local adaptation [2]. Depending on its relationship to the environment, phenotypic variation may be either adaptive or non-adaptive. Adaptive phenotypic variation often occurs between populations that live in different environments and is associated with local adaptation [1,3]. Phenotype-environment correlations have been documented in a wide variety of taxa from plant [4] to invertebrates [5,6] and vertebrates [7][8][9], particularly with respect to the morphology-environment correlation. Functionally important morphological traits that are highly associated with reproductive success, heat exchange, water transfer and locomotion are particularly suitable to studies of speciation and population evolution [10].
Integrative analyses that combine molecular phylogeny, phylogenetic biogeography and phenotypic evolution represent a powerful approach to identify divergent clades with or without phenotypic differentiation, to detect population genetic structure, and to assess early stages of the speciation process [5,11,12]. Studies on lizards have showed that use of different habitats may lead to divergent selection on traits that define body size, body shape, coloration pattern and/or scale characteristics (size, number and scutellation), resulting in morphological diversification among populations or species [13][14][15]. However, to date, few studies have used an integrative approach to address morphological and species diversification of lizards in the Old World.
Toad-headed lizards of the reproductively bimodal genus Phrynocephalus (Agamidae) inhabit desert, arid and semiarid regions in Central and West Asia and North-Northwest China, with all viviparous species restricted to the Qinghai-Tibetan Plateau (QTP) and adjacent regions ( Fig. 1) [16]. The eastern part of the range of the viviparous taxa is occupied by a group of three closely related but taxonomically controversial species, P. guinanensis, P. putjatia and P. vlangalii [17][18][19][20]. Phrynocephalus vlangalii is the most widespread species and inhabits arid and semiarid habitats in the western part of the group's range across an altitudinal range from 2200 to 4500 m, P. putjatia is the oldest species restricted to steppe desert habitats at relatively low altitudes (2200-3300 m) around Qinghai Lake, and P. guinanensis is the most narrowly distributed species restricted to sand dunes (2700-3500 m) in the south of Qinghai Lake ( Fig. 1) [17,[19][20][21][22]. The ranges of P. putjatia and P. vlangalii overlap around Qinghai Lake, but neither occurs in the range of P. guinanensis [17,20]. Morphological data support a valid species of P. guinanensis, but genetic data do not support that distinction [17,19]. So the currently accepted status of P. guinanensis is an ecotype of P. putjatia [19].
However, as the ecotype hypothesis has yet to be empirically tested, a knowledge gap remains. In order to fill the gap we collected specimens from 28 localities (Fig. 1), we downloaded climatic data form WorldClim and trimmed to each sampling locality, took morphological measurements and used molecular markers (mitochondrial ND4 gene and nine microsatellite loci) to assess the structure and clustering of specimens. We then calculated distances based on all of them and compared the dissimilarity matrices. We aim to explore possible scenarios that may explain the discordance between genetic and morphological patterns, and to test whether morphological divergence among these species is associated with local adaptation. We predict that, if the ecotype hypothesis were true, the morphology should be well correlated with climate, or at least more than with genetics.

Genetic polymorphism
We obtained a sequence of 684 base pairs (bp) of the mitochondrial ND4 gene, which contained nine singleton variable sites and 133 parsimony informative sites.  Table 1 for detailed information on sampling localities Eight haplotypes were shared by two species (four by PG and PP, two by PP and PV, and two by PG and PV), and only one haplotype was shared by all three species (Additional file 5: Figure S1). Within individual localities, haplotype diversity varied from 0 to 0.82, and nucleotide diversity from 0 to 0.052 (Table 1). Haplotype diversity (h ± SD) was 0.92 ± 0.01 in PG, 0.94 ± 0.01 in PP, and 0.86 ± 0.02 in PV; nucleotide diversity [(π ± SD) × 10 3 ] was 7.48 ± 4.04 in PG, 23.07 ± 11.44 in PP and 64.23 ± 31.19 in PV (Table 1).
A total of 462 lizards were genotyped and scored at nine microsatellite loci. The number of alleles per locus varied from 14 to 60, with a mean of 37. The mean observed heterozygosity was 0.582, and the mean expected heterozygosity was 0.916 (Additional file 2: Table S2).

Relationships among mtDNA haplotypes
A clade of PV included individuals from PV16, PV20 and PV29; individuals from PV17 formed a clade; individuals from PV18 were admixed in branch. Because of low support values at several nodes, there was a polytomy consisting of PP, PV17, and all other groups of the admixed section, within which a clade of PP included individuals from seven localities (PP22-27) northeast of Qinghai Lake, PV17 and the remaining haplotypes of the three species. Approximately 89% (25/28) of PG haplotypes formed a PG subclade (Fig. 2). Median-joining network based on ND4 haplotypes showed similar grouping patterns to the gene tree and all clades and subclades were recovered (Additional file 5: Figure S1). The mean pairwise distance was 2.0% between PG and PP, 7.1% between PG and PV, and 7.3% between PP and PV. The mean pairwise distance within species was 0.8% for PG, 2.4% for PP, and 7.1% for PV.

Population structure
Assignment tests based on nine microsatellite loci identified two distinct genetic clusters (Fig. 3a). One (red) groups individuals from all three species together, and the other (green) groups individuals from PP and PV. Two major genetic clusters were revealed in PP with one including individuals from localities northeast of Qinghai Lake, and the other including individuals mostly from localities south of the lake. Individuals of PG showed a pure genetic cluster, while individuals of PV had admixed assignment (Fig. 3b). At larger values (3-4) of K, additional clusters appeared. When STRUCTURE was run under the assumption that the data represented three separate populations (K = 3), individuals from localities south of Qinghai Lake were still assigned to their respective clusters (green), but individuals from localities south of the lake and PG individuals were assigned to two groups with moderate probability (red and blue), PV individuals were assigned to a distinct, third cluster (blue) (Fig. 3b).

Morphological divergence
All examined morphological variables except tail length differed among the three species (Table 2). All 11 mensural variables differed between the sexes; only two (superciliaris and dorsal scales) of 11 meristic variables differed between the sexes ( Table 2). PCAs on the body dimensions and scalation characters performed separately for each sex showed that mean scores on the first two axes differed among the three species in both sexes (Additional file 3: Table S3, Additional file 6: Figure S2). Mean scores on the first axis differed among the three species, and in both sexes mean scores were greatest in PG and smallest in PV (Table 3, Fig. 4). Differences were also found between mean scores on the second axis in females, with mean scores being greatest in PV and smallest in PG (Table 3, Fig. 4).

Climatic differences
A PCA of nine climatic variables for 28 localities revealed that the first two components accounted for 80% of the variance (Additional file 4: Table S4). Mean PC scores on the first axis (F 2, 25 = 10.20, P < 0.001; PG a , PP b , PV b ) differed significantly among the three species, while mean PC scores on the second axis did not (F 2, 25 = 1.29, P = 0.29). Overall, climatic differences were more evident between PP and PG than between any other pairs of species (Fig. 5).

Distance correlation analysis
The first morphology PC axis (M1) was positively related to the first climate PC axis (C1) in males (F 1, 19 = 18.22, P < 0.001), and so was in females (F 1, 19 = 18.24, P < 0.001) (Fig. 6). The single Mantel tests for the combined data matrix showed that: (1) geographic distance was significantly related to C1 and M1 in both sexes, and to genetic distance inferred from the ND4 gene; and (2) the first climate PC axis was significantly related to the first morphology PC axis in both sexes (Table 4). In both sexes M1 was significantly related to genetic distance inferred from the ND4 gene. Morphological divergence and genetic distance were spatially patterned and both were climatically dependent (Table 4). Holding C1 constant with the partial Mantel test, we found that the coefficient of a correlation between morphological divergence and genetic distance was 0.187 for males and 0.176 for females, and in both sexes the correlation was not statistically significant (Table 4). Holding geographic distance

111-173
Data are analyzed using two-way ANOVA (SVL and meristic variables) or two-way ANCOVA (mensural variables other than SVL) with SVL as the covariate and species and sex as the factors PG P. guinanensis, PP P. putjatia, PV P. vlangalii, F females, and M males constant, we found once again that C1 significantly correlated with M1 in both sexes (Table 4).

Discussion
Lineage separation and divergence form a temporal process in which populations may accumulate genetic, ecological, and/or morphological changes which make organisms better adapt to their environments, until eventually they are reproductively isolated and form separate species [8,23]. Viviparous Phrynocephalus species form a monophyletic lineage that diverged from the oviparous taxa 9.78 Ma, with the most recent common ancestor of viviparous species dated to 5.04 Ma [24]. The species studied herein do not occur in syntopy (Fig. 1). Of these species, only PP has a range overlapping with oviparous congenerics, largely because it evolved earlier than did other viviparous Phrynocephalus species currently found on the QTP and at relatively low altitudes allowing oviparous reproduction [22,24,25]. The divergence of these species from other viviparous Phrynocephalus species on the QTP is dated to 3.79 ± 0.67 Ma, while the earliest speciation event within the complex is dated to 3.09 ± 0.61 Ma [22,24], following the recent uplift of the QTP (3.6-0.01 Ma). It is therefore likely that environmental changes accompanied by the uplift of the QTP, imposed strong selective forces on local Phrynocephalus populations, and promoted morphological and species diversification. In this study, we found weak genetic differentiation but pronounced morphological divergence between species, and that the morphological diagnoses of species boundaries were not supported by genetic evidence. From this study we can draw the following conclusions. First, PG, PP and PV are not reciprocally monophyletic (Fig. 2). Second, morphological divergence is climatically (ecologically) rather than genetically dependent (Fig. 6). Third, PG is genetically and morphologically more similar to PP than to PV (Figs. 2, 3 and 4).

Weak genetic divergence
Genetic divergence inferred from the ND4 gene was correlated not only with the first climate PC axis (C1) but also with geographic distance (Table 4). This finding allows us to conclude that geographic distance and environmental humidity (or aridity) have major roles in driving genetic divergence between species. In the mtDNA tree, although each species is wildly polyphyletic, we can see similar genetic distance corresponding to the groups identified by morphological characters (Fig. 2). Results of the single Mantel test show a significant correlation between morphological divergence and genetic divergence inferred from the ND4, and significant climatic correlates of morphological and genetic divergence (Table 4). It is worth noting, however, that the morphological-genetic correlation disappeared when holding C1 constant (Table 4). This finding, together with the result that M1 was significantly correlated with C1 in both sexes, indicates that climatic (ecological) dissimilarity rather than genetic divergence has a key role  Microsatellite-based population genetic analyses showed considerable population level admixture. Twenty-nine out of 462 individuals could be assigned to one of the two identified groups with lower than 70% probability, which supports the occurrence of historical introgressive hybridization at the nuclear genetic level. The unclear assignment between PP individuals from south of Qinghai Lake and PG might result from a lack of geographical barriers, fast and recent population expansion, relatively homogeneous habitat, or a combined effect of these factors.
We found two main monophyletic mtDNA clades that separate populations of PV16, PV20 and PV29 to the rest populations. Noble and his colleagues [20] also found two deeply diverged clades in both mtDNA and nuclear markers that were largely congruent with PV and PP; however, there are many individuals with a nuclear genome composition from one species while with mtDNA haplotype from another in ten sampling sites. The admixed mtDNA clade and the individuals sharing the same haplotypes between species confirm the occurrence of historical introgressive hybridization events between species [20]. Two major genetic clusters in PP were found, respectively corresponding to the Qinghai Lake Basin and the southeast of this basin [19]. PG is genetically very close to PP, as revealed by the fact that the mean pairwise distance between PG and PP was only 2.0% (Figs. 2 and 3). In the mtDNA tree, the lack of resolution of star-shaped clade suggests that these groups diverged quite rapidly. Low genetic diversity and clear pure genetic clustering suggest that PG divergence was a very recent event, presumably as a consequence of adapting to desert environments resulting from the uplift of the plateau. High haplotype diversity and low nucleotide diversity indicate rapid recent population expansion in PG. Additionally, the PCA of climatic variables revealed significant climatic niche separation between species (Fig. 5). This result supports the idea that spatially heterogeneous natural selection can lead to population divergence and ecological speciation even in the presence of gene flow [23]. The admixed mtDNA clade and the individuals sharing the same haplotypes between species imply the occurrence of historical introgressive hybridization events between species. Many hybrids (29/462) with an admixed nuclear genome were detected, and we can expect the presence of individuals with admixed or hybrid genomes as a consequence of hybridization events. In addition, high levels of gene flow between three species suggest that these species may suffer from hybridization. Taken together, the three lines of genetic analyses (mtDNA, STRUCTURE and microsatellite based estimations of migration) all suggest ongoing gene flow between species.

Adaptive morphological evolution
Species inhabiting different habitats may experience phenotypic divergence in a suite of traits as a result of adaptation to divergent environments [26]. Using different habitats may lead to divergent selection on a number of fitness-related morphological traits, and the morphologyenvironment correlation has been identified in a number of lizard species [13,14]. For instance, on the gypsum sand dunes of White Sands, data across three different lizard species show that morphological traits are under strong and multifarious selection, and present evidence of the essential factors for divergence [27].
In this study, morphological differences are evident and show adaptive divergence in response to local environments (habitat type in particular). Similar to the pattern of variation in scale number or size reported for lizards of the genera Anolis [28] and Sceloporus [14], our data show that species in more arid environments have fewer larger (inferred from the inverse relationship between scale size and number) scales to reduce skin exposure and thus the amount of evaporative water loss ( Table 2). Scale number is a heritable trait that is likely to respond to ecologicallybased natural selection pressures along environmental gradients, with the complexity of scale hinges, the surface area of skin, and thus the capacity of heat and water exchange increasing with scale number [14,27,29]. In agreement with earlier studies of the species studied herein [17,21], our data show that these species differ morphologically from each other, with body size and scale characters accounting best for morphological divergence between species. Of the three morphological groups, PG and PV are most completely separated, with PP in between (Table 3, Fig. 4). Morphologically, all specimens could be clearly assigned to the species recently described [17,21]. However, morphological divergence is incongruent with genetic divergence inferred from both mitochondrial and microsatellite DNA data sets, as revealed by the fact that the three species do not form any clear lineages or genetic clusters that can be assigned to individual species already described [17,21].
Climatic PC scores differed among the three species. Holding geographic distance constant using the partial Mantel test, we found a significant correlation between ecological divergence (climate PC1) and morphological divergence (morphology PC1) ( Table 4). These findings suggest that morphological differences between species result from local adaptation. Different habitats can generate strong divergent selection and allow adaptive divergence in space even if gene flow is initially substantial [9,[11][12][13]. It is therefore likely that these species exhibit morphological divergence due to their differences in habitat preference. Morphological divergence could restrict gene flow, such as by sexual selection linked to morphological traits or coloration [23,30]. Initial restriction on gene flow could enhance further divergence, and then generate reproductive barriers [31]. Ecological divergence acts as a legitimate isolating mechanism reducing the rate of recombination between divergent habitat types [2], and can therefore drive the evolution of additional intrinsic isolating mechanisms through reinforcement [30].

Species differentiation process
Environmental changes on the QTP, especially desertification and landcover change, have likely driven population divergence and promoted the speciation of the previously so called P. vlangalii (including PP) at 2.29 Ma. After the uplift of the QTP about 1.7 Ma, many areas on the plateau rapidly became arid and some lakes began to disappear [32]. Qinghai Lake became very large during the Middle Pleistocene, caused by the violent lift of the plateau, and then remained stable for a long period due to a drying climate in the Late Pleistocene and declination of its water level since the Holocene [33]. The spreading deserts might have forced P. guinanensis to adapt to sand dunes colonized the adjacent area, resulting in the secondary contact of these three species and potential hybridization. Subsequently, recent gene flow may result in convergence at neutral loci, whereas divergent ecology and selection maintain adaptive differences in morphology.

Conclusions
Our data show that body size and scale characters account best for morphological divergence between species in this group of Phrynocephalus lizards. Morphological divergence is related to habitat types that differ climatically. Morphologically, the species boundary is clear between P. guinanensis and P. vlangalii but not between other two pairs of species. Weak genetic differentiation and pronounced morphological divergence could have resulted from high levels of gene flow and historical introgressive hybridization between species that live in different environments. Our study supports the idea that natural selection in spatially heterogeneous environments can lead to population divergence even in the presence of gene flow. Our study provides a better understanding of genetic, morphological and ecological divergence among closely related species using different habitats, and reveals the initial adaptation to different environments.

Animal collection and treatment
All procedures described in this study were approved by the Animal Care and Use Committee of Nanjing Normal University (2011-04-008). We collected lizards between May and July 2011 from 28 localities around Qinghai Lake (Fig. 1, Table 1). Nine of these localities are occupied only by P. guinanensis (PG), 14 by P. putjatia (PP), and five by P. vlangalii (PV). We identified species based on diagnostic characters reported for these three species [17,21,34]. A total of 175 PG (89 females and 86 males), 195 PP (95 females and 100 males) and 106 PV (54 females and 52 males) adults were used for the collection of morphological (11 mensural and 11 meristic characters) data (Table 2). Morphological information for each species × sex × sampling locality combination with a sample size ≥5 was provided in Additional file 1: Table S1. Following the collection of morphological data, the most distal 2-3 mm of the tail tip was excised from each lizard. Lizards were then released at their site of capture. Tissue samples preserved in absolute ethanol were deposited at Nanjing Normal University.

Mitochondrial DNA amplification and sequencing
Total genomic DNA was extracted from each of 426 individuals using standard phenol-chloroform methods [35]. A partial sequence of the mitochondrial ND4 gene was amplified using forward (ND4) and reverse (Leu) primers [36]. Thermal cycling was performed with initial denaturation for 5 min at 95°C followed by 35 cycles for 50 s at 95°C, 45 s at 58°C, 1 min at 72°C and a final extension for 10 min at 72°C. PCR products were purified and sequenced with each of the PCR primers on an ABI 377 sequencer.

Genetic polymorphism
For mitochondrial DNA data we calculated the number of segregating sites, haplotype diversity and nucleotide diversity for each population (locality) and all populations combined using DnaSP 5.10.1 [41]. Fu's (1997) Fs [39] and Tajima's (1989) D [40] were used to detect departures from the mutation-drift equilibrium that could indicate past demographic changes or selection.
For microsatellite DNA data, parameters such as the number of alleles per locus, average allelic richness, observed heterozygosity, expected heterozygosity, the Hardy-Weinberg equilibrium, and exact tests of linkage disequilibrium between pairs of loci for each population were calculated using ARLEQUIN 3.5 [42] and FSTAT 2.9.3.2 [43].

Phylogeography and population structure
Sequences were aligned using Clustal_X 1.81 [44] with default parameters, and then optimized by eye in MEGA 5 [45]. Mean sequence divergences among major clades were calculated using MEGA 5 and the pairwise Kimura two-parameter (K2P). Bayesian phylogenetic analyses were performed using MrBayes 3.1.2 [46]. We used two oviparous Phrynocephalus species as outgroups: P. albolineatus (GenBank Accession No. AY054002) and P. axillaris (HM235646). Three partitions (the three codon positions of the ND4 sequence; 1st: HKY + G, 2nd: HKY + G and 3rd: GTR + I) were applied to the data and models of molecular evolution were selected for each partition using MrModeltest 2.3 [47]. Four Markov Chains Monte Carlo (MCMC) chains were run for 2.0 × 10 7 generations. Two independent runs were performed to allow additional confirmation of the convergence of MCMC runs. Two runs from random starting trees resulted in the same topology with negligible differences in clade credibility values. We used NETWORK 4.6.1.0 [48] to generate a median-joining network for all individuals of the three species. To facilitate data presentation and interpretation, we used an initial star-contraction procedure with a star connection limit of two to reduce the data set [49].
We examined each population's demographic changes by calculating the raggedness index of the observed mismatch distribution according to the population expansion model implemented in ARLEQUIN 3.5 [42]. We used parametric bootstrapping (1000 replicates) in ARLEQUIN 3.5 to test the goodness-of-fit of the observed mismatch distribution. Whether regional or pooled samples matched the spatial expansion model was estimated by the sum of squared deviations statistic.
We used STRUCTURE 2.3.2 [50] to identify genetically distinct groups among microsatellite genotypes with a burn-in of 5 × 10 7 and 5 × 10 8 iterations without prior population information, following the admixture model. We conducted 10 replicate runs for each specified value of K (the most likely number of populations) from 1 to 20. Individual assignment probability, LnP(D) and convergence between runs were used to assess the most likely value of K, and the most likely number of clusters was estimated according to Evanno and his colleagues [51].

Morphological analyses
We used two-way ANOVA [for snout-vent length (SVL) and all meristic variables] or ANCOVA (for other mensural variables with SVL as the covariate) to examine morphological differences between sexes and among species. Prior to parametric analyses, data were tested for the homogeneity of variances using the Bartlett's test, and for the normality of data using the Kolmogorov-Smirnov test. Tukey's post hoc test was performed on the traits that differed among species. We performed a principal components analysis (PCA) on 22 morphological variables to show positions of three groups (PG, PP and PV; see Fig. 1 for abbreviation definitions) on a two-dimension plane. For the variables that were related to body size, we removed the size effect by using residuals from the regressions against SVL. All statistical analyses were performed with Statistica 10.0 (Tulsa, OK, USA). Throughout this paper, descriptive statistics are presented as mean ± SE and range, and the significance level is set at P < 0.05.

Ecological divergence
In order to evaluate ecological distinctiveness among groups, we used ArcGIS 10.1 to extract values of the 19 climatic variables available in the WorldClim database (http://worldclim.org/version2) at 30 arc-seconds resolution [52]. In order to remove the effect of colinearity, we performed pairwise correlation comparisons between 19 bioclimatic variables and used nine variables that were not highly correlated (r < 0.85) in subsequent analyses. We performed a PCA on the nine climatic variables to reduce the number of predictor variables in our data set, and plotted PC scores on a two-dimension plane. We used linear regression analysis to test if morphology PC scores were correlated with climate PC scores.

Distance correlation tests
We performed a series of single and partial Mantel tests using the 'vegan' package 2.3-5 in R [53] with significance determined using 1000 permutations, testing the correlation between various dissimilarity matrices: (1) geographic distance, (2) climatic PC scores, (3) morphological PC scores, and (4) genetic distance based on mtDNA (genetic distance 1) and microsatellites (genetic distance 2). The climatic and morphological PC scores were converted into distance matrices in R. Only samples that had data for all variables were included in single and partial Mantel tests. Euclidian distance matrices (matrix 2 and 3) were compiled in Statistica 10.0 using the clade means calculated from PC scores for each lizard on a two-dimension plane. Pairwise Fst values [54] were estimated as the matrix of genetic distance using both mitochondrial and microsatellite data, with the procedure implemented in ARLEQUIN 3.5 [42]. Specifically, we tested the following two hypotheses after accounting for geographic distance between localities: (1) climate predicts morphological differences; and (2) climate predicts genetic divergence.