The origin of multiple clones in the parthenogenetic lizard species Darevskia rostombekowi

The all-female Caucasian rock lizard Darevskia rostombekowi and other unisexual species of this genus reproduce normally via true parthenogenesis. Typically, diploid parthenogenetic reptiles exhibit some amount of clonal diversity. However, allozyme data from D. rostombekowi have suggested that this species consists of a single clone. Herein, we test this hypothesis by evaluating variation at three variable microsatellite loci for 42 specimens of D. rostombekowi from four populations in Armenia. Analyses based on single nucleotide polymorphisms of each locus reveal five genotypes or presumptive clones in this species. All individuals are heterozygous at the loci. The major clone occurs in 24 individuals and involves three populations. Four rare clones involve one or several individuals from one or two populations. Most variation owes to parent-specific single nucleotide polymorphisms, which occur as heterozygotes. This result fails to reject the hypothesis of a single hybridization founder event that resulted in the initial formation of one major clone. The other clones appear to have originated via post-formation microsatellite mutations of the major clone.


Introduction
Questions concerning origin and evolution of parthenogenesis in lizards have received considerable attention in recent years [1][2][3][4][5][6]. Parthenogenetic lizard species are ideal organisms for studying the genetic and ecological basis of hybridogeneous speciation, understanding the mechanisms of hybrid parthenogenesis, generation and evolution of genetic diversity, determination of possible factors of natural selection, and hybrid disfunction caused by genetic incompatibilities of admixed genomes or re-established balance of cytonuclear genomes. Unisexual species have been found in less than 0.1% of all vertebrate species [2].
Darevskia rostombekowi (spelling following Murphy [35]) is one of the seven parthenogenetic species that arose from the interspecific hybridization of D. portschinskii (paternal species) and D. raddei (maternal species). The parental species belong to different clades of Darevskia [13,24] and the formation of parthenoforms appears to be phylogenetically constrained [14]. Phylogenetic factors also restrict the formation of some other parthenogenetic species [4]. Darevskia rostombekowi have a chromosome set of 2n = 38 [36], are characterized by fixed heterozygosity of allozyme loci [34] and exhibit low variability of restriction sites of mitochondrial DNA inherited from D. raddei [24]. Studies of 35 allozyme loci from populations in northwestern and central Armenia revealed no allozyme variability, which suggested that this species could be monoclonal [30]. This was unlike the parthenogenetic species D. dahli, D. armeniaca, and D. unisexualis, in which several allozyme clones were found [31]. At present, the paternal species D. portschinskii occurs in southern Georgia, northern Armenia and western Azerbaijan, and the maternal species D. raddei is distributed in Armenia and adjoining Georgia, western Azerbaijan, and northern/western Iran. Several discrete populations of D. rostombekowi occur in Lori and Tavush provinces of northern Armenia (elevation range 700-1500 m), and one small, isolated high-montane (2000 m) population along the southeastern shore of Lake Sevan (Gegharkunik province), in the vicinity of the village of Zagalu (renamed recently to Tsovak). The Sevan and Vardenis mountain ranges, which rise to more than 3000 m, isolate this population from the nearest populations inhabiting the foothills near the Kura River by Lake Gai-Gel', Azerbaijan. The mountains are insurmountable geographical barriers for D. rostombekowi. No allozyme analysis exists for this population.
Analysis of mitochondrial DNA (mtDNA) showed that all northern matrilines of Armenian D. rostombekowi (Papanino, Gosh, Spitak) originated from one southern Armenian population (Yeghegnadzor) of D. raddei [24]. Initial analyses of mtDNA (cytochrome b; Cytb) from D. rostombekowi at Tsovak indicated a maternal origin from D. raddei at Yeghegnadzor [37], but subsequent study based on other mtDNA and nuclear DNA markers [38][39][40] resolved a single origin. All specimens from Tsovak had a single nucleotide substitution in Cytb, which suggested the existence of a second matriline of D. rostombekowi specific to Tsovak. In addition, interpopulation comparisons of genetic variation in D. rostombekowi, based on multilocus DNA fingerprinting with microsatellite hybridization probes, and RAPD markers for the nuclear genome, indicated significant differences between lizards from northern Armenian and Tsovak [41,42]. These findings, plus considerable morphological variation observed among lizards from different populations [13], highlight new questions concerning intra-and inter-population genetic variation and clonal diversity.
Our study examines the clonal diversity in D. rostombekowi including variation at three microsatellite-containing loci identifies genotypes of D. rostombekowi from four Armenian populations. We apply a new approach to genotyping the related parthenogenetic species D. dahli [43] and show that D. rostombekowi-previously presumed to be a monoclonal parthenogenetic species-has several clones. New, the approach reveals parent-specific markers consisting of microsatellites and single nucleotide polymorphisms (SNPs) located out of the microsatellite cluster for each allele of each locus. The SNPs yield direct information about interspecific hybridization founder events, and microsatellite variability provides information concerning possible mutations in the initial hybrid clones that give rise to new genotypes. Interspecies allele comparisons involve homologous microsatellite loci from the parental bisexual lizard species D. raddei and D. portschinskii. All alleles present in D. rostombekowi also occur in the parental species.

Materials and methods
DNA samples of D. rostombekowi (n = 42) from four populations, D. raddei (n = 65) from 13 populations, and D. portschinskii (n = 27) from two populations in Armenia were analyzed (Table 1). Sample localities were shown in Fig 1. All DNA was extracted from samples obtained about 15 years ago. Parthenogenetic and bisexual lizards were collected between 1997 and 2006 from their natural habitats in Armenia. All collecting sites occurred on governmental land that did not involve protected areas. Therefore, specific permissions were not required for collecting in these locations. Our field studies predated the listing of D. rostombekowi as an Endangered Species by the IUCN (2009) and the species is 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. All specimens were captured with a noose. 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.
Three loci, Du215, Du281, and Du323 were PCR-amplified using previously described primer pairs [43,44]. These loci were chosen because each contained a tetranucleotide microsatellite cluster that provided sufficient resolution of the individual allelic PCR products in non-denaturing polyacrylamide gel electrophoresis to allow the direct sequencing of the individual alleles. In general, the procedure for the isolation and sequencing of individual alleles was carried out as described previously [44,45].
Data on genetic variation of these loci in D. rostombekowi [46] and on allelic variation of the homologous loci in D. raddei and D. portschinskii were shown in S1 Table. 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, and 48˚C for Du323); and 40s at 72˚C followed by one cycle of 5 min at 72˚C. PCR products (15 μl) were loaded onto an 8% non-denatured polyacrylamide gel (to separate allelic variants for each locus) and run for 12 h at 60V. A 50bp 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 (KM573728-KM573762; HM014002-HM014003; KR559279-KR559316). The number of alleles, allelic richness, 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 Fstat v.2.9.3.2 [47], GenePop v.4.2, and Web-version of POPTREEW [48]. Analysis of molecular variance AMOVA was carried out using the package poppr [49]. We implemented multiple alignments of microsatellite allele sequences obtained from the three species by using MUSCLE [50]. These matrices were used for the building a neighbor joining network with MEGA v.6.0.6 [51]. The method of coding genotypes was described previously [45].

Results
All individuals of D. rostombekowi were heterozygous at the three loci and contained two alleles that differed from each other in length and structure of the microsatellite clusters. Two of them, Du281 and Du323, had also single nucleotide variants (SNVs) in fixed positions of the flanking regions (S1 Table). Unlike SNPs in sexual species, SNVs in hybrid genomes of parthenogenetic lizards result from the codominant loci of divergent parental genomes. They represent parent-specific markers in parthenogenetic species. Thus, distinct clonal SNVs would likely have owed to independent hybridization founder events. The number of alleles in D. rostombekowi varied from two to five depending on the locus (S1 Table). The allelic variants of Du281 and Du 323 formed distinct groups according to parent-specific SNVs. SNVs in alleles 1-3 of Du281 formed the SNV set TATA, whereas those in allele 4 of Du281 formed the set CGCG. Similarly, SNVs in alleles 1 and 2 of Du323 had two sets of SNVs: CT and AC. These variants reflected the hybrid origins of the parthenogenetic species and the heterozygosity of its genome. Du215 did not possess parent-specific SNVs in the vicinity of the microsatellite cluster (S1 Table). The dinucleotide insertion TT (-98/99) Du215(rost)3 occurred only in one individual from Gosh and this insertion was not found in those populations of the parental species analyzed (see below and S1 Table). In D. rostombekowi (n = 42), Du215 had five alleles, Du281 had four, and Du323 had two (S1 Table).
The origin of SNVs at Du281 was shown schematically in S1 Fig based on alignment of the microsatellite sequences and neighboring regions. The SNVs of Du281 connected with parental alleles at the heterozygous locus.
To identify genotypic diversity, we constructed allelic combinations in the 42 individuals of D. rostombekowi. Based on both single nucleotide and microsatellite variants, analyses resolved five genotypes that differed in their frequencies and geographic distribution ( Table 2).
Individuals with identical genotypes were assumed to represent distinct clonal lineages. Occurring in 24 individuals (57.1% of the total cohort), genotype 1 was most abundant; it was found in three populations. Genotypes 2 and 3 were represented by eight individuals each (19% of the total cohort) in two and one population, respectively. Finally, rare genotypes 4 and 5 were represented by one individual each (2.38% of the total cohort). The genotypic diversity among the four populations of D. rostombekowi varied from 0% to 75% ( Table 2). The highest level of genotypic diversity was observed at Gosh, where rare genotypes 4 and 5 occurred; both variants might have originated via post-formation microsatellite mutations from the most abundant, and presumptively initial, genotype 1. The lowest level of genotypic diversity was observed at Tsovak, which had genotype 3 only. This genotype differentiated remote Tsovak from all populations in northern Armenia. We carried out AMOVA of four D. rostombekowi populations. However, the results of AMOVA were insignificant probably because of low sample size in most of the localities. Population genetic indices for k loci and genotypes 1-5 were given in Table 3. Expected heterozygosity varied from 0.500 to 0.750 (average, 0.534) and the number of alleles ranged from two to four (average, 2.7). Values of allelic richness ranged from 1.995 to 4.000 (average, 2.620). The highest allelic scores, number of alleles and allelic richness occurred at Gosh for the all three loci.
Indices of gene diversity in males and females of the parental species D. raddei (n = 65) and D. portschinskii (n = 27) were given in Table 4. Their allelic diversity was higher than that of the parthenogenetic populations. Du323 had lowest number of alleles, ranging from 2 to 7, depending on species. The number of alleles for Du281 varied from 8 to 15, and for Du215 from 9 to 17, depending on species. The values of allelic richness at all three loci were equal for males and females in D. raddei, and a little more for male D. portschinskii. Du215 and Du281 had heterozygote deficits. That of Du281 also occurred in male D. raddei. Female D. raddei had an excess of heterozygotes at Du323. Deviations from Hardy-Weinberg expectations probably owed to small sample sizes. The heterozygote deficits may have owed to a Wahlund effect due to the 65 individuals of D. raddei coming from 13 localities.
To determine if D. rostombekowi originated from a single or multiple interspecies hybridization event(s), genotype-specific markers, formed by combinations of parent-specific SNVs of both alleles in Du281 and Du323 were evaluated (Fig 2; S1 Table). All genotypes matched the parent-specific combinations TATA/CGCG (Du281) and CT/AC (Du323). This suggested a common origin of all five genotypes owing to a single interspecies hybridization event. However, we could not rule out the possibility of independent crossings of parental individuals possessing the same alleles in these loci. Genotype 5 in one individual from Gosh had the same genotype-specific markers for loci Du281 and Du323, but differed in having unique allele Du215(rost3) with its dinucleotide insertion TT(-98/99). This insertion was not detected among analyzed parental alleles (see below; S1 Table) and, thus, it probably represented a de novo mutation. Genotypes 1-5 differed from each other by microsatellite sequences only. Some of the rarer genotypes may have arisen through post-formation microsatellite mutations of the dominant genotype. Genotype 3 from Tsovak was the most divergent. Genotypes 4 and 5 appeared to be independent events and they were observed only in the Gosh population, which had the most variability of genotypes. No unique genotypes were found at Papanino and Spitak. The analysis of the spatial-frequency distributions of the more widely distributed genotypes 1 and 2, and population-specific genotype 3 revealed a dependence between the frequencies of the clones and the geographical disjunction between the three northern populations and the Tsovak population.

Discussion
Given that lizards with identical genotypes represent a distinct clone, the several microsatellite clones of D. rostombekowi reject the hypothesis of monoclonality based on allozymes [30]. However, it corresponds with the considerable amount of morphological variation in this species [13]. The observed difference in patterns of allozymes and microsatellites likely reflects the less variable nature of allozymes compared with microsatellites. Parthenogenetic D. armeniaca, D. dahli and D. unisexualis have more than one allozyme clone [52] and microsatellites may reveal additional clones. Our approach, which is based on detection of both, microsatellite variability and single nucleotide variation of DNA sequences outside of microsatellite clusters, promises to reveal even higher levels of clonal diversity in these species. For example, an analysis of 35 allozyme loci revealed five clones in D. dahli, whereas three microsatellite-containing loci detected 11 clones [43]. Identical genotype-specific markers occur all clones of D. rostombekowi, suggesting that it was formed by one hybridization event. Identical SNVs but different microsatellite markers show that the hybridization event involved very similar female D. raddei and male D. portschinskii in one population. It remains unclear which of these clones of D. rostombekowi formed first. Assuming that unisexual species originating from a single hybridization exhibit a widespread and common clone with a few rare clones [23], genotype 1 of D. rostombekowi might be ancestral ( Table 2). The other clones likely represent post-formation microsatellite mutations. The highly unstable (GATA) n -containing locus in parthenogenetic D. unisexualis appears to have accumulated mutant alleles in some offspring of the first generation; de novo mutations occur via deletion or insertion of a single microsatellite repeat [45]. Such mutations might be responsible for the origin of genotypes 4 and 5 of D. rostombekowi at Gosh and genotype 2 at Papanino and Spitak (Table 2). Surprisingly, major genotype 1 and genotypes 2, 4, and 5 remain undetected at Tsovak. This discovery questions the origin of genotype 3. Genotype 3 may owe to an independent hybridization event. However, we do not have a sufficient number of genotype-specific markers to test for this possibility. Alternatively, genotype 3 may have originated from common genotype 1, only to completely replace it. Ecological, physical and climatic differences among the habitats of various populations of D. rostombekowi may probably explain interpopulational differences in the distribution of various microsatellite genotypes [53][54][55]. Diversification may also relate to the frequency of interspecific mating at different localities, the relative fitness of clonal lineages and/or outcompeting hybrids its sexual progenitors [56].
Analysis of mitochondrial DNA shows that D. raddei from Yeghegnadzor is the matriarch of all northern populations (Papanino, Gosh, Spitak). However, the finding of a second mitotype of D. rostombekowi specific for the Tsovak population [37] indicates the possible participation of different females of D. raddei in the hybridization. To examine this possibility, additional genotype-specific markers from more microsatellite loci are necessary.
Parthenogenetic species of Darevskia exhibit low levels of clonal diversity and a high degree of mtDNA similarity with respect to their parental species. Their recent origins appear to best explain these patterns, while they may have largely different ages. Darevsky [57] suggested that all parthenogenetic species of Darevskia originated in the Caucasus after the last glacial maximum and that all species have had the same amount of time to evolve new clones, approximately 5000-7000 years ago. However, a multilocus analysis for the parthenogens D. unisexualis, D. bendimahiensis and D. uzzeli dated their origins to between 200 000 and 70 000 years ago [6]. These results further support the independent origins of parthenogenesis in Darevskia and these predate the last glacial maximum. Different parental lineages hybridized several times in different geographic regions [14]. Further, Irwin et al. [58] provided an estimate of a divergence rate of approximately 10% per million years for the silent substitution at the third-codon position, based on mammalian Cytb data. If this calculation is applicable to lizards, the age of D. rostombekowi would be approximately 200 000 years old (2% divergence) [59].
Electrophoretic studies of proteins and mitochondrial DNA have shown that parthenogenetic lizards tend to be multiclonal. A high level of variability of allozyme loci was discovered in Heteronotia binoei (Gekkonidae) [16]. Genetic heterogeneity was also demonstrated for the American parthenogenetic teiid lizards Aspidoscelis tessellatus and A. neomexicanus [60] and in parthenogenetic D. dahli, D. armeniaca, and D. unisexualis [31]. Darevskia rostombekowi has two times lower clonal diversity than D. dahli [43]. Although data are limited, this suggests that either D. rostombekowi may have a more recent origin than D. dahli, or it suffered a severe bottleneck.
Some cytogenetic and molecular characteristics of the genome of hybrid species differ from those of the bisexual species due to disturbance in the segregation of chromosomes. This results in cytogenetic variability, different genomic activity and mutations in some loci [14,31,42,61]. Karyotypic instability of hybrid genomes may lead to the genetic diversity of parthenogenetic species and the appearance of chromosomal mutations [36]. Such mutations may govern the formation of new clones and/or geographically isolated chromosomal races. Our analyses cannot rule out the possibility that the isolated population of D. rostombekowi inhabiting the coast of Lake Sevan constitutes a chromosomal race and/or a clone that emerged because of chromosomal mutations occurring in the course of karyological evolution.
In summary, our analyses discover multiclonal genetic structure in the parthenogenetic lizard D. rostombekowi. Future studies may determine whether the clonal diversity owes to a single or multiple interspecific hybridization events, and identify the role mutations play in altering the initial clone(s).