- Research article
- Open Access
Multiple interspecific hybridization and microsatellite mutations provide clonal diversity in the parthenogenetic rock lizard Darevskia armeniaca
BMC Genomics volume 19, Article number: 979 (2018)
The parthenogenetic Caucasian rock lizard Darevskia armeniaca, like most other parthenogenetic vertebrate species, originated through interspecific hybridization between the closely related sexual Darevskia mixta and Darevskia valentini. Darevskia armeniaca was shown to consist of one widespread allozyme clone and a few rare ones, but notwithstanding the origin of clonal diversity remains unclear. We conduct genomic analysis of D. armeniaca and its parental sexual species using microsatellite and SNP markers to identify the origin of parthenogenetic clonal lineages.
Four microsatellite-containing loci were genotyped for 111 specimens of D. armeniaca, 17 D. valentini, and four D. mixta. For these species, a total of 47 alleles were isolated and sequenced. Analysis of the data revealed 13 genotypes or presumptive clones in parthenogenetic D. armeniaca, including one widespread clone, two apparently geographically restricted clones, and ten rare clones. Comparisons of genotype-specific markers in D. armeniaca with those of its parental species revealed three founder-events including a common and two rare clones. All other clones appeared to have originated via post-formation microsatellite mutations in the course of evolutionary history of D. armeniaca.
Our new approach to microsatellite genotyping reveals allele-specific microsatellite and SNP markers for each locus studied. Interspecies comparison of these markers identifies alleles inherited by parthenospecies from parental species, and provides new information on origin and evolution of clonal diversity in D. armeniaca. SNP analyses reveal at least three interspecific origins of D. armeniaca, and microsatellite mutations in these initial clones give rise to new clones. Thus, we first establish multiple origins of D. armeniaca. Our study identifies the most effective molecular markers for elucidating the origins of clonal diversity in other unisexual species that arose via interspecific hybridization.
Naturally occurring unisexual reproduction occurs in less than 0.1% of all vertebrate species [1, 2] and true parthenogenesis has been described only in squamates [1, 3,4,5]. Recently, the origin and evolution of parthenogenesis, in particular in lizards, has received considerable attention [6,7,8,9,10,11,12,13]. Unisexuality usually arises through interspecific hybridization between genetically related sexual species [14, 15] and introgressive hybridization is common among species of lizards Jancuchova-Laskova et al. . Regardless, hybridization between the presumed ancestral species of parthenogenetic lineages of Caucasian rock lizards (Darevskia) does not necessary result in parthenogenesis . Among Squamata, many parthenogenetic species are triploid (e.g., some Aspidoscelis), however, all unisexual species of Darevskia (Lacertidae) are diploid. Due to their hybrid origin, these species have the genetic diversity of their parental sexual progenitors and commonly exhibit fixed heterozygosity at codominant loci . Sister-chromatid pairing maintains heterozygosity, which is critical for offsetting the reduced fitness in parthenogenetic species .
Genetic studies of unisexual vertebrate species rely on assessments of intraspecific variability and clonal diversity. Their genetic diversity may owe to having the original clones but from different founders, post-formation mutations (especially in hypervariable microsatellite loci), and genetic recombination in rare cases of subsequent crossings [19, 20]. Range-size, the age of species and some environment factors may also affect clonal diversity [21,22,23].
Darevsky (1958) first discovered parthenogenesis in the family Lacertidae . Later, other parthenogenetic species were described [25, 26]. Darevskia has at present 24 sexual and seven parthenogenetic species. [11, 14, 25]. The phylogeny of Darevskia reconstructed using mitochondrial DNA sequences and allozyme data  includes the D. caucasica, D. saxicola and D. rudis major groups of sexual species. The formation of parthenospecies is constrained phylogenetically [7, 27]. They arose from hybridization of females in the D. caucasica group and males in the D. rudis group. Analysis of DNA fingerprint markers has identified genetic polymorphism in all parthenogenetic species of Darevskia [28,29,30,31,32].
Darevskia armeniaca has a wide distribution involving central Armenia, southern Georgia, northeastern Anatolia and northwestern Azerbaijan [33,34,35]. They live at elevations between 800 and 2700 m [33, 36]. This species arose from the hybridization of Darevskia mixta (D. caucasica group) and Darevskia valentini (D. rudis group) [14, 15, 23, 27, 37, 38]. Darevskia mixta, which occurs in the eastern part of the Meskheti Range, Lesser Caucasus Mountains and on southern slops of the Greater Caucasus between the valleys of the Rioni and Khobi rivers, is endemic to Georgia. The central part of the Lesser Caucasus was the most likely area of origin for D. mixta . In comparison, D. valentini occurs in eastern Anatolia, Armenia, and adjoining Georgia [15, 33] at elevations between 1900 and 3110 m. In Armenia, this species is locally abundant in montane habitats . Triploid hybrids between D. armeniaca and sexual D. valentini occur in some Armenian populations [34, 40, 41]. The zone of sympatry in Kuchak has an extremely high number of hybrids . The majority of them are sterile triploid females, but males with fully developed reproductive systems and presumably fertile females also exist . Putative fertile triploids probably exist in D. unisexualis, and in some other parthenogenetic species of reptiles [41,42,43,44]. In most cases, spatial gaps separate the ranges of sexual and unisexual populations of Darevskia . Unisexual forms have a superior competitive ability to their parental forms in a given spatio-temporal setting , expansion of the parthenogens could cause retreat of parental species from the contact zones, preventing further hybridization events .
Recently, Tarkhnishvili et al. (2017) postulated alternative hypothesis of origin of D. armeniaca. Apparently, it arose from backcrosses of male D. valentini with parthenogenetic D. dahli . Mitochondrial DNA analysis confirmed data of Fu et al. (1999) that D. armeniaca and D. dahli descend maternally from D. mixta from a limited geographic area in central Georgia . The majority of both parthenogens shared the same genotypes at two microsatellite loci, but they differed at the other three loci used. Therefore, they suggested that the origin of D. armeniaca first included hybridization between D. mixta and D. portschinskii followed by backcrosses of these parthenogens with male D. valentini . However, new studies are needed to determine which of two or both scenarios of D. armeniaca origin have matter.
To test the alternative hypotheses, we examine the clonal diversity and its origin in D. armeniaca using 111 samples and microsatellite genotyping and single nucleotide polymorphisms (SNPs) located outside of each microsatellite cluster [46, 47]. The SNPs data yield direct information about interspecific hybridization founder events, and microsatellite variability provides information about possible mutations in the initial clones. Further, analyses use partial sequences from the mitochondrial gene encoding cytochrome b (MT-CYTB) from 14 Armenian populations of D. armeniaca plus the parental species.
DNA samples of D. armeniaca (n = 111) from 14 populations, D. valentini (n = 17) from four populations in Armenia, and D. mixta (n = 4) from one population in Georgia were analyzed (Table 1). Sample localities were shown in Fig. 1.
All DNAs were isolated from blood samples of lizards. Lizards were captured by Danielyan F.D. and Arakelyan M.S. 15–20 years ago with a noose. These species are not protected by CITES. The study was approved by the Ethics Committee of Moscow State University (Permit Number: 24–01) and was carried out in strict accordance with their ethical principles and scientific standards. Blood samples were taken from the tail veins of lizards under chloroform anesthesia, and then these lizards were released. DNA was isolated from lizard blood by using the standard phenol–chloroform extraction method with proteinase K, and resuspended in TE buffer, pH 8.0.
Four loci, Du215, Du281, Du323, and Du47G were PCR-amplified using previously described primer pairs [46, 48, 49]. Data on genetic variation of these loci in D. armeniaca and on allelic variation of the homologous loci in D. valentini and D. mixta  were shown in Additional file 1: Table S1. PCR was performed on 50 ng of DNA in a total volume of 20 μl using a GenePak PCR Core Kit (Isogene) and 1 μM of each primer. The reaction conditions were as follows: one cycle of 3 min at 94 °C; 30 cycles of 1 min at 94 °C; 40s at the annealing temperature (58 °C for Du215, 50 °C for Du281, 52 °C for Du323, and 54 °C for Du47G); and 40s at 72 °C followed by one cycle of 5 min at 72 °C. PCR products (7 μl) were loaded onto an 8% non-denatured polyacrylamide gel (to separate allelic variants for each locus) and run for 12 h at 60 V. A 100 bp ladder (Fermentas) was used as a size marker. The amplified products were visualized by staining DNA in the gel with ethidium bromide. Well-resolved individual PCR products, which corresponded to the two individual alleles of the locus, were excised from the gel, purified by ethanol precipitation, and sequenced directly in both directions using a chain termination reaction with an ABI PRISM BigDye Terminator v.3.1 on an Applied Biosystems 3730 DNA analyzer. Allelic identity was checked and confirmed via the comparison of sequences obtained independently. All unique de novo sequences were deposited in GenBank (GU972533-GU972535; HM070259-HM070264; HM013992-HM013994; KT070998-KT071004; GU972551; GU972553; HM013997; KM573717-KM573727; MH187988-MH187999). Sometimes for genotypes 9–13, which were represented by only one individual, the PCR products of the Du47G locus were excised from the polyacrylamide gel, purified and cloned into pMos blue vectors according to standard procedures (pMos blueBlunt ended Cloning kit RPN 5110, Amersham Biosciences). The clones were amplified in MOSBlue competent cells grown at 37 °C, and sequenced as described previously .
Diversity parameters for each locus were estimated only for D. armeniaca, because sample sizes of the parental sexual species were insufficient. The number of alleles, allelic richness (Rs), as a measure of allele counts adjusted for sample size, and expected heterozygosity, as a measure of gene diversity, were calculated per locus and per population by using R package Poppr v.2.5.0 , GenePop v.4.2, and Web-version of POPTREEW .
A statistical parsimony haplotype network was used to visualize variation. It was calculated using TCS software v.1.21, with gaps being considered as a second state . This approach has been used at the population-level for comparing mitochondrial SNPs, which have linear arrangements . Clonality in parthenogenetic lizards resulted in homologous alleles having linear arrangements with little or no recombination. This served to link the parthenogenetic genotypes with mutations (repeat number changes). Tarkhnishvili et al. (2017) used a similar approach for linking microsatellite genotypes of D. dahli and D. armeniaca from Georgia, Armenia and Turkey, but with Network v.5.0 [13, 55]. The method of coding genotypes was described previously .
A 320 bp fragment of MT-CYTB was amplified and sequenced for 30 specimens of D. armeniaca (2–3 individuals from each population) including those with a distinct microsatellite genotype, and four specimens of D. mixta from one population. The primers used were L14841: 5′-CCATCCAACATCTCAGCATGATGAAA-3′ and H15149: 5′-GCCCCTCAGAATGATATTTGTCCTCA-3′ , as described for these taxa in . PCR analysis was performed as described in . The amplicons were sequenced on an automated sequencer (Applied Biosystems 3730 DNA analyzer). Sequence alignment was performed with BioEdit v.7.0 .
All 111 specimens of D. armeniaca, including 10 individuals with coloration difference  from Artavaz (Hankavan), 17 specimens of D. valentini, and four D. mixta were analyzed successfully using locus-specific PCR and DNA sequencing of PCR amplificants. The 320 bp fragment of MT-CYTB did not vary among D. armeniaca and D. mixta; all sequences assigned to haplotype A of D. mixta .
All but four individuals of D. armeniaca (from Dsegh, Dilijan-Semyonovka Pass and Stepanavan) were heterozygous at microsatellite loci used; the alleles differed from each other in length and structure, and in single nucleotide variations (SNVs) in fixed positions of the flanking allelic regions (Additional file 1: Table S1), which was similar to previous reports for D. dahli, D. unisexualis, and D. rostombekowi [32, 46,47,48]. Locus Du323, also had a unique (AC)n microsatellite cluster that differentiated parental alleles.
In D. armeniaca, Du215 and Du323 had three alleles, Du281 had four, and Du47G had seven (Additional file 1: Table S1). Allelic variants formed distinct groups according to the fixed SNVs, which in hybrid genomes resulted from the different parental genomes. They corresponded to parent-specific markers such that unique clonal SNVs likely reflected independent founder events. SNVs in alleles 1 and 2 of Du215 formed the set TGC and allele 3 of Du215 formed the set ACT. Similarly, SNVs in alleles of Du281 had T and C variants, and SNVs in alleles of Du323 had CT and AC variants. The seven alleles of Du47G formed three sets of SNVs: TAGT, TTCA, and AAGA. Each sets associated with specifically organized microsatellite clusters. Thus, together the sets of SNVs and microsatellite clusters differentiated distinct genotypes inherited by parthenospecies from parental species.
Analysis of the parental species showed homozygosity for Du215 and Du323 in D. mixta (n = 4), and Du281 and Du47G had two alleles each. In D. valentini (n = 17), Du215 was homozygous, Du281 had five, Du323 had six, and Du47G had 10 alleles (Additional file 1: Table S1). All parental alleles contained microsatellite clusters and variable nucleotides (SNVs) at fixed positions in neighboring regions which differentiated D. mixta and D. valentini.
Microsatellite motifs and/or the sets of SNVs discerned parental alleles. Most alleles of Du215, Du281, Du323, and Du47G in D. armeniaca occurred in the parents. Alleles Du215(arm)1 and 2 coincided with allele Du215(mix)1, and allele Du215(arm)3 coincided with Du215(val)1. Further, alleles Du281(arm)1 and 2 coincided with Du281(mix)1 and 2, and alleles Du281(arm)3 and 4 coincided with Du281(val)1–5. Alleles Du323(arm)1 and 2 with the specific (AC)6 microsatellite cluster occurred among Du323 alleles from D. valentini, and allele Du323(arm)3 with specific (AC)5 microsatellite cluster coincided with Du323(mix)1. Finally alleles, Du47G(arm)3–6 corresponded to alleles of Du47G in D. valentini, and Du47G(arm)7 coincided with Du47G(mix)2. Due to the high mutation rate of microsatellite DNAs , in some cases we did not find direct correlation between distinct alleles of D. armeniaca and its parents. Parental alleles with the SNV TAGT, which was specific for Du47G(arm)1 and 2, was not found probably due to either genome divergence through genetic recombination or mutation in D. armeniaca, or insufficient parental sampling. Allelic combinations of four loci identified genotypic diversity in D. armeniaca (Fig. 2).
Analyses resolved 13 genotypes that differed in their frequencies and distribution (Table 2). Widespread genotype 1 occurred in 61 individuals (54.9% of the total cohort). Genotype 4 (n = 10; 9%) and genotype 6 (n = 5; 4.5%) were found in four and three populations, respectively. All other genotypes (2, 3, 5, 7–13) occurred in one or two populations (n = 35; 32.5%). Although relatively uncommon, genotypes 2–5 made up each the majority of individuals in some populations: genotype 2 at Harich (12 of 18 specimens), genotype 3 at Kuchak (5 of 8), genotype 4 at Pushkin Pass (5 of 7), and genotype 5 at Lchashen (7 of 9) (Table 2). All 10 color-variant individuals from Artavaz (Hankavan) had genotype 1.
Genotypic diversity of D. armeniaca varied from 0 to 66.7% (Table 2). The highest level of genotypic diversity was observed at Dsegh and Sotk, which had two genotypes in three individuals. Five genotypes, one common and four rare, occurred in the nine individuals from Stepanavan. The lowest levels of genotypic diversity were observed at Sevan, Lchashen, Artavaz (Hankavan), and Tezh (Pambak Ridge), which had genotype 1 only.
Population genetic indices for four loci and genotypes 1–9 were given in Table 3. Five individuals with unique genotypes 10–13, and populations at Sevan and Lchashen, which were represented by one individual each, were excluded from the analysis. The estimates of expected heterozygosity varied from 0.51 to 0.67 (average 0.55–0.63 depending on locus) whereas observed heterozygosity did not vary among loci and populations (Table 3). The number of alleles varied from 2 to 5 (average, 2.00–2.83 depending on locus). Values of allelic richness ranged from 1.89 to 3.20 (average, 1.94–2.24 depending on locus). The highest values of heterozygosity occurred at Dsegh (loci Du215, Du281, Du323) and Sotk (locus Du47G). The highest values of allelic richness occurred at Dsegh (loci Du215 and Du281), Sotk (loci Du215, Du281, and Du47G) and Medved-gora (vicinity of Stepanavan) (locus Du323), while the highest values of allelic number varied from 2 to 5 depending on the locus in the populations studied.
Population genetic indices for four loci and four populations (17 individuals) were given in the Additional file 2: Table S2. The estimates of expected heterozygosity varied from 0.50 to 0.89 (average 0.52–0.84 depending on locus), whereas observed heterozygosity from 1 for locus Du215 (all populations) to 0.33–0.94 for other loci. The number of alleles varied from 2 to 6 (average, 2–5 depending on locus). Values of allelic richness ranged from 1.79 to 3.43 (average, 1.96–3.36 depending on locus). The highest value of allelic richness occurred in the population at Lchashen (loci Du47G, Du323, Du215). The highest values of expected heterozygosity occurred in Adis (locus Du47G). These data demonstrated that higher genetic variability in populations of D. valentini rather than populations of D. armeniaca. Linkage disequilibrium analysis (Additional file 3: Table S3) suggested that all populations except Lshachen had no association between loci. The population at Lshachen appeared to have linked loci probably because of either a recent bottleneck or association of microsatellite markers with unknown loci under selection.
Combinations of parent-specific SNVs and those of D. armeniaca differentiated between single or multiple interspecies hybridization event(s). The structural composition of the 13 genotypes (Additional file 1: Table S1) were shown schematically in Fig. 2. Genotypes 1–9 matched all parent-specific SNV combinations TAGT/TTCA (Du47G), TGC/ACT (Du215), T/C (Du281), and CT/AC (Du323). This did not reject the hypothesis of a common origin from a single hybridization event. However, the analysis did not rule out the possibility of independent crossings of the parental individuals because they differed from each other by microsatellite sequences only at loci Du47G, Du281, and Du323. Some of the rarer genotypes may have arisen via post-formation microsatellite mutation.
Genotypes 10–12, which occurred in four individuals from three populations, matched parent-specific SNV combinations at all loci, but they differed from genotypes 1–9 by the parent-specific SNV combination TTCA/TTCA for locus Du47G. These four specimens were uniquely homozygous at Du47G and differed from each other only by microsatellite sequences at Du215 and Du281. Therefore, this rejects the hypothesis of a single origin of D. armeniaca, and their variation likely arose through microsatellite mutations. Further, rare genotype 13 (n = 1) differed from all other genotypes by the SNV combination TTCA/AAGA at locus Du47G. This also rejected the null hypothesis. Consequently, analyses pointed to at least three independent hybridization events in the genesis of the 13 genotypes of D. armeniaca.
The TCS network displayed the geographical distribution of genotypes 1–9 (Additional file 4: Figure S1). The network displayed a star-like appearance that was typical of mutations deriving from a central, ancestral genotype.
In their review of allozyme variation in parthenogenetic lizards, Parker et al. (1989) proposed that species having a single origin will usually have a widespread numerous clone along with a few rare clones, where species with multiple hybrid origins are highly variable with random allele combinations . In some cases, high diversity in allozymes but low variation in mitochondrial DNA suggests multiple origins from a geographically restricted sample of females [22, 58, 59]. Add to this, low diversity in allozymes and mitochondrial DNA suggests restricted origins, both geographically and numerically [60, 61]. In a laboratory creation of hybridogenetic Poeciliopsis lineages, relatively low clonal diversity resulted from selection of the most-fit clones from a broad spectrum of genotypes [2, 16].
The clonal diversity of parthenogenetic species was detailed using allozyme and mitochondrial DNA analysis [16, 23, 27, 38, 45, 62,63,64,65]. Allozyme analyses resolved several clones in all species except for D. rostombekowi), and all consist in one major widespread clone and a few rare clones. Some peculiarities in allozyme clone patterns were found for D. armeniaca [38, 66].
Darevskia armeniaca has a karyotype of 2n = 38  characterized by fixed heterozygosity of allozyme loci [38, 66] and exhibits no  or low variability of mitochondrial DNA inherited from D. mixta . Using five loci, Uzzell and Darevsky (1975) found one clone in 11 specimens of D. armeniaca from two populations from Armenia . MacCulloch et al. (1995) examined 35 loci and 75 specimens from seven Armenian populations and found one widespread clone with two rare clones . One of the rare clones occurred in 19 out of 27 individuals at Dilijan (Papanino) and the other one the only clone in Kuchak (n = 2). Fu et al. (2000) used 35 allozyme loci and 117 specimens, including some lizards with the color variation described by Darevsky (1992) and Danielyan (1999) as being four clones, one being common [37, 66, 68]. One rare allozyme clone dominated in two populations and a rare color-variant clone at Ankavan (n = 2) differed at two loci . This is unlike D. dahli, whose morphological difference did not correspond with either allozyme or microsatellite markers [16, 46].
MacCulloch et al. (1995) attributed the clonal variation that associated with morphological data in D. armeniaca to mutations [37, 38]. They argued that (1) some rare clones had alleles not detected in the parental species; (2) the pattern of clonal variation in D. armeniaca was typical of parthenogenetic lizards of single hybridization origin such as in teiid Aspidocelis ; and (3) overall variation in D. armeniaca was less than that found in species of multiple hybrid origin, such as parthenogenetic gecko Heteronatia binoei . Alternatively, Fu et al. (2000) suggested that although mutation was a possible explanation of the origin of the clonal variation in D. armeniaca, but multiple origins was equally likely . The rare clones in the Kuchak and Dilijan (Papanino) may have owed to an independent origins because of dominance of rare alleles and (2) the peripheral distribution to the common clone suggested multiple origin . Further, the young age of D. armeniaca implied from low mtDNA variation and substantial allozyme variation favored the multiple-origin scenario over a rapid accumulation of mutations .
Our analyses indicate that at least three interspecific hybridizations were involved in the genesis of D. armeniaca and microsatellite mutations followed these. A previous allozyme analysis reported four clones with the rare ones occurring in restricted populations . This allozyme pattern generally followed Parker et al.’s (1989) model , which is consistent with single hybridization origin. Nevertheless, Fu et al. (2000) noted that multiple origins could also explain the clonal diversity .
Our modified microsatellite genotyping  detects both microsatellite and SNVs variability. The method reveals 13 genotypes in 111 individuals using four microsatellite loci. Previously, the approach revealed a high level of clonal diversity (11 clones) in D. dahli  and rejected the hypothesis of monoclonality of D. rostombekowi , and our analyses also detect greater variation in D. armeniaca than did allozymes. Thus, this method is a more precise measure of assessing genetic variability compared with allozymes.
Previously, the two color-variant individuals of D. armeniaca from Ankavan were investigated using allozymes . These specimens differed from the common allozyme clone forming separate clones. However, our data failed to distinguish 10 specimens of D. armeniaca with differences in coloration; they all have the most abundant and widespread genotype 1. The color-variants consist of about half of the samples from Artavaz (Hankavan). Similarly, our microsatellite genotyping data , as well as allozyme  and mtDNA  studies, failed to distinguish the color-varieties of parthenogenetic D. dahli.
The pattern of distribution of the clones rejects the hypothesis of a single origin for D. armeniaca. All individuals with genotypes 1–9 exhibit an identical combination of parent-specific SNVs, suggesting their common origin via one hybridization event. Given that D. mixta and D. valentini are the parental species for D. armeniaca [14, 23, 37, 38, 66], the identical SNVs but different microsatellite markers indicates that the hybridization event involved one population. The five individuals with genotypes 10–12 also exhibit an identical combination of parent-specific SNVs, suggesting their common origin, but one of that differs from individuals with genotypes 1–9. Clonal diversity in both groups owes to microsatellite unstable (GATA)n mutations only. Such mutations in parthenogenetic D. unisexualis occur in one generation via deletion or insertion of a single repeat at one or at both alleles of the locus . Finally, the individual from Harich has genotype 13 and a unique combination of parent-specific SNVs at locus Du47G. Allele Du47G(arm)7 was inherited by D. armeniaca from D. mixta (Additional file 1:Table S1). Therefore, this individual may represent a third hybridization event.
The rather high clonal diversity in D armeniaca likely owes to its multiple hybrid origin and subsequent microsatellite mutations. Similarly, most of the observed clonal diversity within parthenogenetic D. dahli in 9 out of 11 detected genotypes owe to microsatellite mutations within the common clone, and two out of 111 individuals were suggested to be members of independent hybridization events . Analyses of parthenogenetic D. rostombekowi  revealed one common and four rare clones, each represented by one or several individuals from one or two populations. The results were consistent with single hybridization origin of D. rostombekowi, with clonal diversity arising via post-formation microsatellite mutations .
Tarkhnishvili et al. (2017) hypothesized that D. armeniaca originated from a series of crosses between D. valentini and parthenogenetic D. dahli, rather than between D. valentini and D. mixta . Their microsatellite genotyping of Du215, Du281, Du481, Du323, and Du47G revealed identical genotypes for D. dahli and D. armeniaca at Du323 and Du47, but not at the other three loci. They suggested that backcrosses and mutations best explained the genetic and clonal diversity than a multiclonal origin. However, species of Darevskia exhibit rather high levels of genomic similarity and very similar microsatellite alleles. They assumed alleles of equal length had identical nucleotide sequences. Given our discoveries, more detailed molecular data on the structure of the alleles and genotypes is needed to support their hypothesis. Our discovery of multiple origins provides a more parsimonious hypothesis than the complex hypothesis of backcrossing. Further, projections on the historical distribution of sexual and parthenogenetic species during glacial waves is necessary to document the possibility of hybrids between D. dahli and D. valentini, as well as possible hybridizations between D. mixta and D. portschinskii that might be according to  the first stage in origin of D. armeniaca.
Parker et al. (1989) proposed that unisexual species originating through a single hybridization event will exhibit a common clone with a few rare clones . Accordingly, among clones 1–9, genotype 1 might be ancestral (Table 2). The TCS network (Additional file 4: Figure S1) has a star-like structure with common genotype 1 occupying the central location and the others differing by from one to two mutational events. This star-like structure is consistent with a recent origin and diversification of clones. The same inference is not possible for genotypes 10–12 because none of them is numerous and widespread. Dsegh and Harich do not appear to have genotype 1. However, genotype 11 at Dsegh and genotype 13 at Harich arose independently from genotype 1. Genotype 4 at Dsegh and Harich, and genotypes 2, 6 and 7 at Harich appear to be variants of genotype 1. If so, then their occurrence in these localities must owe to the dispersal of lizards from other regions.
MacCulloch et al. (1995) and Fu et al. (2000) observed some peculiarities in allozyme-clone variation [62, 65]. Among their five populations of D. armeniaca, rare clones made up the majority of the individuals in Kuchak and Dilijan (Papanino), although the common clone covered most of the species’ distribution. This contributed to the suggestion that numerous rare clones might indicate independent hybrid origins . However, our parent-specific SNVs for these clones are identical to those of genotype 1.
The distribution of microsatellite clones in D. armeniaca is not exceptional. Among 11 detected clones in the parthenospecies D. dahli , two were numerous and widespread geographically and all others were rare, although one rare clone in “Dendropark” (near Stepanavan) occurred in 6 out of 9 individuals. Similar analysis revealed one common and four rare clones in D. rostombekowi, and one of the rare clones occurred in all 8 individuals at Tsovak . Thus, in all three parthenogenetic species investigated, the numerous rare clones most likely arose via post-formation microsatellite mutations of the common clone, and not through independent interspecific hybridizations. Ecological differences among the habitats of various populations, and the relative fitness of clonal lineages of D. armeniaca and other parthenogenetic Darevskia may explain the differences in the distribution and existence of successful rare clones.
In summary, the analyses suggest that the clonal diversity in D. armeniaca derives from three interspecific hybridizations and subsequent microsatellite mutations. Future studies may identify the role mutations play in altering of the initial clones. The methodological approach, which is based on the detection of parent-specific microsatellite and SNV markers, can elucidate the origin of genetic and clonal diversity in other unisexual species that arose via interspecific hybridization.
Our interspecific genomic analysis used microsatellites and single nucleotide polymorphisms (SNP). A comparison of these markers between parthenogenetic D. armeniaca and its parental sexual species D. valentini and D. mixta reveals 13 genotypes or presumptive clones. Clonal diversity in D. armeniaca appears to result from at least three independent interspecific hybridization events. All other clones of D. armeniaca appear to have derived from subsequent microsatellite mutations. This methodological approach, which is based on the detection of parent-specific microsatellite and SNP markers, can be applied for study clonal diversity in other unisexual species that arose via interspecific hybridization.
Vrijenhoek RC, Dawley RM, Cole CJ, Bogart JP. A list of known unisexual vertebrates. In: Dawley RM, Bogart JP, editors. Evolution and ecology of unisexual vertebrates. Albany: New York state museum bulletin 466; 1989. p. 19–23.
Neaves WB, Baumann P. Unisexual reproduction among vertebrates. Trends Genet. 2011;27:81–8.
Avise JC. Clonality: the genetics, ecology and evolution of sexual abstinence in vertebrate animals. Oxford: Oxford University Press; 2008.
Hedges SB, Marion AB, Lipp KM, Marin J, Vidal N. A taxonomic framework for typhlopid snakes from the Caribbean and other regions (Reptilia, Squamata). Caribb Herpetol. 2014;49:1–61.
Booth W, Schuett GW. The emerging phylogenetic pattern of parthenogenesis in snakes. Biol J Linn Soc Lond. 2016;118:172–86.
Manriquez-Moran NL, Cruz FR, Murphy RW. Genetic variation and origin of parthenogenesis in the Aspidoscelis cozumela complex: evidence from mitochondrial genes. Zool Sci. 2014;31:14–9.
Jancuchova-Laskova J, Landova E, Frynta D. Are genetically distinct lizard species able to hybridize? A review. Curr Zool. 2015;61:155–80.
Jancuchova-Laskova J, Landova E, Frynta D. Experimental crossing of two distinct species of leopard geckos, Eublepharis angramainyu and E. macularius: viability, fertility and phenotypic variation of the hybrids. PLoS One. 2015;10(12):e0143630. https://doi.org/10.1371/journal.pone.0143630.
Abdala CS, Baldo D, Juárez RA, Espinoza RE. The first parthenogenetic Pleurodont Iguanian: a new all-female Liolaemus (Squamata: Liolaemidae) from western Argentina. Copeia. 2016;104:487–97.
Deakin J, Edwards MJ, Patel H, O’Meally D, Lian J, Stenhouse R, et al. Anchoring genome sequence to chromosomes of the central bearded dragon (Pogona vitticeps) enables reconstruction of ancestral squamate macrochromosomes and identifies sequence content of the Z chromosome. BMC Genomics. 2016;17:447. https://doi.org/10.1186/s12864-016-2774-3.
Freitas S, Rocha S, Campos J, Ahmadzadeh F, Corti C, Sillero N, et al. Parthenogenesis through the ice ages: a biogeographic analysis of Caucasian rock lizards (genus Darevskia). Mol Phylogenet Evol. 2016;102:112–27.
Freitas S, Vavakou A, Arakelyan M, Dvoretski SV, Crnobrnja-Isailović J, Kidov AA, et al. Cryptic diversity and unexpected evolutionary patterns in the meadow lizard, Darevskia praticola (Eversmann, 1834). Syst Biodivers. 2016;10:184–97.
Tarkhnishvili D, Murtskhvaladze M, Anderson CL. Coincidence of genotypes at two loci in two parthenogenetic rock lizards: how backcrosses might trigger adaptive speciation. Biol J Linn Soc Lond. 2017;121:365–78.
Uzzell ТМ, Darevsky IS. Biochemical evidence for the hybrid origin of the parthenogenetic species of the Lacerta saxicola complex (Sauria, Lacertidae), with a discussion of some ecological and evolutionary implications. Copeia. 1975;2:204–22.
Darevsky IS, Kupriyanova LA, Uzzell T. Parthenogenesis in reptiles. In: Gans C, Billett DF, editors. Biology of Reptilia. New York: John Wiley & Sons; 1985. p. 413–526.
Murphy RW, Darevsky IS, MacCulloch RD, Fu J, Kupriyanova LA, Upton DE, et al. Old age, multiple formations or genetic plasticity? Clonal diversity in the uniparental Caucasian rock lizard, Lacerta dahli. Genetica. 1997;101:125–30.
Simon JC, Delmotte F, Rispe C, Crease T. Phylogenetic relationships between parthenogens and their sexual relatives: the possible routes to parthenogenesis in animals. Biol J Linn Soc Lond. 2003;79:151–63.
Lutes AA, Neaves WB, Baumann DP, Wiegraebe W, Baumann P. Sister chromosome pairing maintains heterozygosity in parthenogenetic lizards. Nature. 2010;464:283–6.
Moritz C, Donnelan S, Adams M, Baverstock PR. The origin and evolution of parthenogenesis in Heteronotia binoei (Gekkonidae): extensive genotypic diversity among parthenogens. Evolution. 1989;43:994–1003.
Cole CJ, Dessauer HC, Barrowclough GF. Hybrid origin of a unisexual species of whiptail lizard, Cnemidophorus neomexicanus, in western North America: new evidence and a review. Am Mus Novit. 1988;2905:1–38.
Dessauer HC, Cole CJ. Diversity between and within nominal forms of unisexual teiid lizards. In: Dawley RM, Bogart JP, editors. Evolution and ecology of unisexual vertebrates. Albany: New York state museum bulletin 466; 1989. p. 49–71.
Parker ED, Walker JM, Paulissen MA. Clonal diversity in Cnemidophorus: ecological and morphological consequences. In: Dawley RM, Bogart JP, editors. Evolution and ecology of unisexual vertebrates. Albany: New York state museum bulletin 466; 1989. p. 72–6.
Moritz C, Uzzel Т, Spolsky C, Hotz H, Darevsky IS, Kupriyanova LA, et al. The maternal ancestry and approximate age of parthenogenetic species of Caucasian rock lizards (Lacerta: Lacertidae). Genetica. 1992;87:53–62.
Darevsky IS. Natural parthenogenesis in certain subspecies of rock lizards (Lacerta saxicola Eversmann). Dokl Akad Nauk SSSR Biol Sci. 1958;122:730–2.
Uetz P. The Reptile Database. (2018). http://www.reptile-database.org. Accessed 02 Oct 2018.
Trifonov VA, Paoletti A, Barucchi VC, Kalinina T, O’Brien PC, Ferguson-Smith MA, et al. Comparative chromosome painting and NOR distribution suggest a complex hybrid origin of triploid Lepidodactylus lugubris (Gekkonidae). PLoS One. 2015;10(7):e0132380. https://doi.org/10.1371/journal.pone.0132380.
Murphy RW, Fu J, MacCulloch RD, Darevsky IS, Kupriyanova LA. A fine line between sex and unisexuality: the phylogenetic constraints on parthenogenesis in lacertid lizards. Zool J Linnean Soc. 2000;130:527–49.
Kan NG, Petrosyan VG, Martirosyan IA, Ryskov AP, Darevsky IS, Danielyan FD, et al. Genomic polymorphism of mini- and microsatellite loci of the parthenogenic Lacerta dahli revealed by DNA fingerprinting. Mol Biol. 1998;32:672–−8.
Tokarskaya ON, Kan NG, Petrosyan VG, Martirosyan IA, Grechko VV, Danielyan FD, et al. Genetic variation in parthenogenetic Caucasian rock lizards of genus Lacerta (L. dahli, L. armeniaca, L. unisexualis) analyzed by DNA fingerprinting. Mol Gen Genet. 2001;265:812–9.
Tokarskaya ON, Martirosyan IA, Badaeva TN, Malysheva DN, Korchagin VI, Darevsky IS, et al. Instability of (GATA)n microsatellite loci in the parthenogenetic Caucasian rock lizard Darevskia unisexualis (Lacertidae). Mol Gen Genomics. 2003;270:509–13.
Malysheva DN, Tokarskaya ON, Petrosyan VG, Danielyan FD, Darevsky IS, Ryskov AP. Genomic variation in parthenogenetic lizard Darevskia armeniaca: evidence from DNA fingerprinting data. J Heredity. 2007;98:173–8.
Ryskov AP. Genetically unstable microsatellite-containing loci and genome diversity in clonally reproduced unisexual vertebrates. Int Rev Cell Mol Biol. 2008;270:319–49.
Arakelyan MS, Danielyan FD, Corti C, Sindaco R, Leviton AE. Herpetofauna of Armenia and Nagoro-Karabakh. Salt Lake City: Society for Study of Amphibians and Reptiles; 2011.
Darevsky IS. Rock lizards of the Caucasus: systematics, ecology and phylogenesis of the polymorphic groups of Caucasian rock lizards of the subgenus Archaeolacerta. Leningrad: Nauka; 1967. [in Russian: English translation published by Indian Natl Sci Doc Ctr, New Delhi, 1978]
Tarkhnishvili D. Evolutionary history, habitats, diversification, and speciation in Caucasian rock lizards. In: Jenkins OP, editor. Advances in zoology research. Hauppauge: Nova Science Publishers; 2012. p. 79–120.
Tarkhnishvilli DN, Gavashelishvili A, Avaliani A, Murtskhvaladze M, Mumladze L. Unisexual rock lizard might be outcompeting its sexual progenitors in the Caucasus. Biol J Linn Soc Lond. 2010;101:447–60.
Darevsky IS. Evolution and ecology of parthenogenesis in reptiles. In: Adler K, editor. Herpetology: current research on the biology of amphibians and reptiles. Proceedings of the first world congress of herpetology. Oxford: Society for the Study of Amphibians and Reptiles; 1992. p. 447–60.
MacCulloch RD, Murphy RW, Kupriyanova LA, Darevsky IS, Danielyan FD. Clonal variability in the parthenogenetic rock lizard, Lacerta armeniaca. Genome. 1995;38:1057–60.
Gabelaia M, Tarkhnishvili D, Murtskhvaladze M. Phylogeography and morphological variation in a narrowly distributed Caucasian rock lizard, Darevskia mixta. Amphibia-Reptilia. 2015;36:4554.
Darevsky IS, Danielyan FD. Diploid and triploid progeny arising from natural mating of parthenogenetic Lacerta armeniaca and L. unisexualis with bisexual L. saxicola valentini. J Herpetol. 1968;2:65–9.
Danielyan F, Arakelyan M, Stepanyan I. Hybrids of Darevskia valentini, D. armeniaca and D. unisexualis from a sympatric population in Armenia. Amphibia-Reptilia. 2008;29:487–504.
Reeder TW, Dessauer HC, Cole CJ. Phylogenetic relationships of whiptail lizards of the genus Cnemidophorus (Squamata, Teiidae): a test of monophyly, reevaluation of karyotypic evolution, and review of hybrid origins. Am Mus Novit. 2002;3365:1–62.
Cole CJ, Hardy LM, Dessauer HC, Taylor HL, Townsend CR. Laboratory hybridization among north American whiptail lizards, including Aspidoscelis inornata arizonae × A. tigris marmorata (Squamata: Teiidae), ancestors of unisexual clones in nature. Am Mus Novit. 2010;3698:1–43.
Wynn AH, Cole CJ, Gardner AL. Apparent triploidy in the unisexual Brahminy blind snake, Ramphotyphlops braminus. Am Mus Novit. 1987;2868:1–7.
Fu J, Murphy RW, Darevsky IS. Limited genetic variation in Lacerta mixta and its parthenogenetic daughter species: evidence from cytochrome b and ATPase 6 gene DNA sequences. Genetica. 1999;105:227–31.
Vergun AA, Martirosyan IA, Semyenova SK, Omelchenko AV, Petrosyan VG, Lazebny OE, et al. Clonal diversity and clone formation in the parthenogenetic Caucasian rock lizard Darevskia dahli. PLoS One. 2015. https://doi.org/10.1371/journal.pone.0091674.
Ryskov AP, Osipov FA, Omelchenko AV, Semyenova SK, Girnyk AE, Korchagin VI, et al. The origin of multiple clones in the parthenogenetic lizard species Darevskia rostombekowi. PLoS One. 2017;12(9):e0185161. https://doi.org/10.1371/journal.pone.0185161.
Korchagin VI, Badaeva TN, Tokarskaya ON, Martirosyan IA, Darevsky IS, Ryskov AP. Molecular characterization of allelic variants of (GATA)n microsatellite loci in parthenogenetic lizards Darevskia unisexualis (Lacertidae). Gene. 2007;392:126–33.
Girnyk АЕ, Vergun АА, Omelchenko AV, Petrosyan VG, Korchagin VI, Ryskov AP. Molecular and genetic characterization of the allelic variants of Du215, Du281, Du323, and Du47G microsatellite loci in parthenogenetic lizard Darevskia armeniaca (Lacertidae). Russ J Genet. 2017;53:472–82.
Badaeva TN, Malysheva DN, Korchagin VI, Ryskov AP. Genetic variation and de novo mutations in the parthenogenetic Caucasian rock lizard Darevskia unisexualis. PLoS One. 2008;3(7):e2730. https://doi.org/10.1371/journal.pone.0002730.
Kamvar ZN, Tabima JF, Grünwald NJ. Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ. 2014;2:e281. https://doi.org/10.7717/peerj.281.
Takezaki N, Nei M, Tamura K. POPTREEW: web version of POPTREE for constructing population trees from allele frequency data and computing some other quantities. Mol Biol Evol. 2014;31:1622–4.
Templeton AR, Crandall KA, Sing CF. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics. 1992;132:619–33.
Posada D, Crandall KA. Intraspecific gene genealogies: trees grafting into networks. Trends Ecol Evol. 2001;16:37–45.
Bandelt HJ, Forster P, Rohl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Kocher TD, Thomas WK, Meyer A, Edwards SV, Paabo S, Villablance FX, et al. Dynamics of mitochondrial DNA evolution in animals, amplification and sequencing with conserved primers. Proc Natl Acad Sci U S A. 1989;86:6196–200.
Hall TA. BioEdit, a user-friendly biological sequence alignment editor and analysis program for windows 95/98/ NT. Nucleic Acids Symp Ser. 1999;41:95–8.
Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9:1657–9.
Densmore LD, Moritz C, Wright JW, Brown WM. Mitochondrial DNA analysis and the origin and relative age of parthenogenetic lizards (genus Cnemidophorus) IV. Nine sexlineatus group unisexualis. Evolution. 1989;43:469–83.
Vyas DK, Moritz C, Peccinini-Seale D, Wright JW, Brown WM. The evolutionary history of parthenogenetic Cnemidophorus lemniscatus (Sauria, Teiidae). II. Maternal origins inferred from mitochondrial DNA analysis. Evolution. 1990;44:922–32.
Moritz C. The origin and evolution of parthenogenesis in Heteronotia binoei: evidence for recent and localized origins of widespread clones. Genetics. 1991;129:211–9.
MacCulloch RD, Fu J, Darevsky IS, Danielyan FD, Murphy RW. Allozyme variation in three closely related species of Caucasian rock lizards (Lacerta). Amphibia-Reptilia. 1995;16:331–40.
MacCulloch RD, Murphy RW, Fu J, Darevsky IS, Danielyan FD. Disjunct habitats as islands: genetic variability in the Caucasian rock lizard Lacerta portschinskii. Genetica. 1997;101:41–5.
Fu J, MacCulloch RD, Murphy RW, Darevsky IS, Kupriyanova LA, Danielyan FD. The parthenogenetic rock lizard Lacerta unisexualis: an example of limited genetic polymorphism. J Mol Evol. 1998;46:127–30.
Fu J, MacCulloch RD, Murphy RW, Darevsky IS, Tuniyev BS. Allozyme variation patterns and multiple hybridization origins: clonal variation among four sibling parthenogenetic Caucasian rock lizard. Genetica. 2000;108:107–12.
Fu J, MacCulloch RD, Murphy RW, Darevsky IS. Clonal variations in the Caucasian rock lizard Lacerta armeniaca and its origin. Amphibia-Reptilia. 2000;21:83–9.
Kupriyanova LA. Genetic diversity of hybrid unisexual species and forms of genus Lacerta (Lacertidae, Reptilia): its possible cytogenetic mechanisms, meiosis of cytogenetics in natural polyploidy forms. Cytologia. 1999;21:1038–46.
Danielyan FD. Natural mutagenesis as a factor of speciation in populations of parthenogenetic lizards. In: Abstracts of the scientific conference devoted to the 75th anniversary of the Department of Zoology of the YSU. Yerevan: Department of Zoology of the YSU; 1998. p. 40–2.
Authors thank to Fedor Osipov for help in preparing the manuscript.
This research was partly funded by RFBR according to the research projects № 17–00-00430 (17–00-00426), № 17–04-00396 and by the RAS Program «Molecular and Cell Biology». The work was conducted on the base of the Center for Collective Use, Institute of Gene Biology, Russian Academy of Sciences (GK-02.451.11.7060). These funding bodies had roles in the experimental work, analysis of the data, and in preparation and writing the manuscript.
Availability of data and materials
All data generated or analyzed during this study are included in this published article and its supplementary information files. All unique de novo sequences were deposited in GenBank (http://www.ncbi.nlm.nih.gov/).
Studied lizards species are not protected by CITES. Lizards were captured by Danielyan F.D. and Arakelyan M.S. 15–20 years ago with a noose. Blood samples were taken from the tail veins of lizards under chloroform anesthesia, and then these lizards were released. The study was approved by the Ethics Committee of Moscow State University (Permit Number: 24–01) and was carried out in strict accordance with their ethical principles and scientific standards.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Allelic variations of microsatellite containing loci in the lizard species D. armeniaca, D. valentini, and D. mixta. (PDF 132 kb)
Table S2. The population indices of gene diversity for four studied loci in twelve sampled populations of D. valentini. (PDF 99 kb)
Table S3. Indices of association in four populations of D. valentini. (PDF 347 kb)
Figure S1. Schematic representation of the TCS network that reflects distribution of genotypes 1–9 in D. armeniaca. Concatenated sequences of D. armeniaca genotypes were analyzed using TCS software v.1.21. Genotypes 10–13 are plotted separately. Population distribution of the genotypes is shown by different colors. Numbers indicate the number of individuals in populations. The black circles show, unsampled, but computer-predicted genotypes. (TIF 802 kb)
About this article
Cite this article
Girnyk, A.E., Vergun, A.A., Semyenova, S.K. et al. Multiple interspecific hybridization and microsatellite mutations provide clonal diversity in the parthenogenetic rock lizard Darevskia armeniaca. BMC Genomics 19, 979 (2018) doi:10.1186/s12864-018-5359-5
- SNP markers
- Clonal diversity