Origin, clonal diversity, and evolution of the parthenogenetic lizard Darevskia unisexualis.

BACKGROUND
The hybridization of female D. raddei and male D. valentini gave rise to the parthenogenetic Caucasian rock lizard Darevskia unisexualis. A previously identified genetic polymorphism in the species consisted of one common and two allozyme clones. Analysis of microsatellites and single nucleotide polymorphisms (SNPs) from the three species yields estimates of clonal diversity and tests the hypothesis of a single origin for D. unisexualis.


RESULTS
Genotyping and sequencing of four microsatellite-containing loci for 109 specimens of D. unisexualis, 17 D. valentini, and 45 D. raddei nairensis identified 12 presumptive clones, including one widespread and 11 rare clones. Most individuals in some localities had a rare clone. Clone-specific alleles in D. unisexualis were compared with those of the parental species. The results inferred a single hybridization event. Post-formation mutations best explain the less common clones.


CONCLUSIONS
Interspecific analyses identify alleles inherited by D. unisexualis from its bisexual ancestors. SNP analyses fail to reject the hypothesis of a single interspecific origin of D. unisexualis, followed by microsatellite mutations in this initial clone. Microsatellites detect higher clonal diversity in D. unisexualis compared to allozymes and identify the likely origins of clones. Our approach may be applicable to other unisexual species whose origins involve interspecific hybridization.

Herein, we use parthenogenetic Darevskia unisexualis to test three hypotheses. 1) The parthenospecies has a single origin from the hybridization of paternal D. valentini and maternal D. raddei [26]. 2) Most clonal variation owes to post-hybridization mutation. And 3) analyses microsatellite loci will detect higher clonal diversity compared to allozymes. To test these hypotheses, assessments of clonal variation are essential. The extent of variation within parthenospecies can depend on the rate of clonal formation [27], ecological specialization of clonal lineages [28], historical biogeography [29], and other processes.
Parthenogenetic D. unisexualis was found to have three allozyme clones, based on analyses of 36 loci from three populations of D. unisexualis in central Armenia (n = 57) [26]. Rare clones occurred in two individuals and all others consisted of a common, widespread clone. Its low level of variation in mitochondrial DNA and allozymes among populations suggested that the founders of D. unisexualis involved very few individuals [25,26]. However, the origin of this variation, whether owing to point mutations, insertion/deletions, multiple origins, or more complex genomic reorganization, remains unresolved, in part due to the species' widespread distribution. The species occurs in East Anatolia and in small, isolated areas in central Armenia (Aragatsotn, Gegharkunik, Kotayk, Lori, and Shirak provinces) [30] where it prefers rocky exposures and its vertical distribution ranges from 1500 to 2300 m a.s.l. Although many populations are large, some are threatened by overgrazing and urbanization. Accordingly, this species is classified as "Near Threatened" by the IUCN and listed as "Vulnerable" in the Red Book of Armenia [30].
Paternal D. valentini occurs in eastern Turkey and high montane habitats (elevations of 1900-3110 m) in central Armenia and adjoining Georgia; populations are locally abundant and IUCN assessed this species as Least Concern [30]. Maternal ancestor D. raddei is widespread throughout central Armenia, with isolated populations in the north and in south-central portions of the country; it also occurs in adjoining Georgia and East Anatolia [31]. Like other congeners, D. raddei prefers stony or rocky habitats at elevations of 1000-2660 m. Individuals are usually abundant and the IUCN assessed it as "Least Concern" [30]. Darevskia raddei has been suggested to be a species-complex, containing the forms "raddei" and "nairensis" whose taxonomic status is still a matter of debate [25,32,33], and this uncertainty extends into the origins of parthenogenetic clones [34]. Notwithstanding, D. raddei nairensis occurs sympatrically with D. unisexualis at Lchap, Armenia (Gegharkunic Province) on the western margin of Lake Sevan [30]. Because the parental species of D. unisexualis exhibit high allozyme variation among populations [35,36], the parthenospecies likely originated from few parental individuals [25]. Analysis of mitochondrial DNA obtained a concordant result; the four populations of D. unisexualis had identical sequences, but populations of D. raddei exhibit variation [25].
Our analyses of D. unisexualis use variation at four microsatellite-containing loci in seven Armenian populations. The same methods were used previously in our assessments of D. dahli [37], D. rostombekowi [38], and D armeniaca [39]. Interspecies comparisons use alleles of homologous loci from D. unisexualis and bisexual parents D. valentini and D. raddei nairensis. Analyses of D. unisexualis and its maternal parent also include partial sequences of mitochondrial cytochrome b (CYTB). Results show that D. unisexualis has a level of clonal diversity similar to ones of other parthenospecies of Darevskia. Analyses provide direct information about interspecific hybridization founder events, and about possible mutations in the initial hybrid clones.

Results
All individuals of D. unisexualis had identical fragments of CYTB. The fragment assigned to haplotype of D. raddei nairensis from Lchashen, Armenia (data not shown; GenBank Accession No. U88613).
Each microsatellite locus in individuals of D. unisexualis had two alleles. Both length and structure of the alleles differed within individuals. Further, the flanking regions of the alleles had single nucleotide polymorphisms (SNPs) in fixed positions ( Fig. 1 and Table S1). All clones had identical combinations of parent-specific SNPs for all loci, which was consistent with an origin from a single interspecies hybridization event; these results did not reject the first hypothesis that a single hybridization event gave rise to D. unisexualis. Further, the alternative hypothesis of multiple origins was rejected because D. unisexualis did not share alleles with multiple variants of either parental species. Assuming this to be true, then we could not reject the second hypothesis that most clonal variation owed to posthybridization mutation.
Variation in SNPs and microsatellites resulted in 12 clones ( Fig. 1 and Table 1). Clone C1 (clone 1) was found in all populations and it occurred in 37 individuals (33.9% total cohort). All other clones occurred in one or two populations only. Clone C2 occurred only at Artavaz (n = 28; 25.7% total cohort), C4 at Noratus only (n = 14; 12.8% total cohort), and C3 dominated at Lchap (n = 14; 12.8% of total cohort). Clones C5-C12 occurred in one or two populations and were found in 1-3 individuals (n = 16; 14.6% total cohort). Clones C10-C12 were represented by only one individual each. Because allozyme analyses only  Table 1 Clones, clone composition, sample size, distribution of clones among populations, and diversity of alleles in D. unisexualis. For clone composition, allelic notation is (allele number in D. nairensis + allele number in D. valentini). Alleles shown in Fig. 1 Clone Clone composition Population Number of individuals (clone frequency) Artavaz Hrazdan Kuchak Lchap Noratus Sevan Tsovak resolved three clones, our microsatellite results failed to reject the third hypothesis that microsatellites will detect more clones than allozymes.
The average values of allelic richness for individual loci varied significantly within the species. Locus Du47G had one allelic variant in all populations of D. raddei nairensis (Table S2), as did Du215 in D. valentini (Table S3). However, the polymorphic loci of the parental species often had greater allelic richness than D. unisexualis. Allelic richness of Du281 in populations of D. raddei nairensis (6.01 ± 0.24) was significantly higher (p < 0.05) than in D. unisexualis (2.42 ± 0.17), as was allelic richness of Du47G (3.36 ± 0.24) in D .valentini versus D. unisexualis (2.10 ± 0.04) (p < 0.05). Nevertheless, the average values of total allelic richness of all loci did not differ significantly among all species (p > 0.05) due to homozygosity of some loci of the parental species.
The TCS network placed the most common clone (C1) in a central location with respect to the other, less common clones. The other clones differed from it by one or two mutations only (Fig. 2). Within populations, clonal diversity in D. unisexualis ranged from 8.0 to 66.7% ( Table 1). The highest levels were observed at Sevan, which had two clones in three individuals, and at Kuchak, in which the 12 individuals had one common and four rare clones. With two clones in 25 individuals, Tsovak had the lowest level of genotypic diversity. The number of alleles varied from 2 to 4 and allelic richness ranged from 1.98 to 3.00 (Table 2). Tsovak (Du281 and Du47G), Hrazdan (Du215 and Du323), and Kuchak (Du323 and Du281) had the highest values of allelic richness.
Population genetic indices for four populations of D. raddei nairensis (42 individuals) were given in Table S2. Populations Ayrivank (n = 2) and Bjni (n = 1) were excluded from analyses owing to small sample sizes. Observed heterozygosity ranged from 0.40 to 1.00 (average 0.60-1.00 depending on locus), and this was similar to expected heterozygosity, which ranged from 0.32 to 0.83 (average 0.32-0.80 depending on locus). From 1 to 10 alleles were observed, depending on locus and population. Depending on locus and population, allelic richness varied from 1 to 6.67. Du281 at Pyunik had highest value of allelic richness. Expected heterozygosity was greatest in Lchashen (Du281). Population genetic indices for D. valentini (17 individuals) were calculated previously [39] and are presented in the Table S3.

Discussion
Within D. unisexualis, the two rare allozyme clones were hypothesized to have resulted via post-formation mutation of the preexisting common clone [26]. This  [24] that little allozyme and mtDNA variation exists in species having a single hybridization origin. Spatially, one widespread clone is common, but a few rare clones also exist. This pattern holds for other parthenogenetic Caucasian rock lizards [40,41], as well as other parthenogenetic lizards [42]. Comparatively, parthenogenetic D. dahli and D. armeniaca also have one common clone and several rare ones [40]. Only parthenogenetic D. rostombekowi exhibited a single allozyme clone [43]. The level of diversity in D. unisexualis was also similar to that found in parthenogenetic Aspidoscelis neomexicanus [44], which also has a hybrid origin. In contrast, parthenogenetic Heteronotia binoei was reported to have much higher variation [42].
Our microsatellites and SNPs reveal higher levels of clonal diversity in parthenospecies of Darevskia than allozymes did. Analyses involving 109 individuals of D. unisexualis from seven populations in Armenia identify 12 clones that differ in their frequencies and population distribution. Analyses of 35 allozyme loci in parthenogenetic D. dahli [40], D. rostombekowi [43], and D. armeniaca [45] resolved five, one, and four clones, respectively, while our genomic approach resolves 11 clones in D. dahli [37], five in D. rostombekowi [38], and 13 in D. armeniaca [39]. Thus, assessments of microsatellites discover more variation than allozymes.
Analyses cannot reject the hypothesis of one hybridization event [26] forming D. unisexualis due to identical SNPs in flanking regions of the microsatellite loci and a single mitochondrial haplotype. However, it remains possible that ancestral parents experienced back crossings. Clones C1-C12 differ from each other only by microsatellite sequences. Variation in lengths of microsatellite alleles surely owes to the high rate of indels, which can occur in one generation [46]. However, future analyses of additional populations of D. unisexualis can test for the possibility of multiple origins, which could account for geographic patterns of alleles, as opposed to mutations within geographic regions. It is particularly important to sample populations from geographically distant locations in Turkey to differentiate between scenarios of dispersal and multiple origins.
Our analyses detect high genotypic diversity in parthenogenetic D. unisexualis similar to those found in parthenogenetic D. dahli and D. armeniaca. Variation in clone frequency could result from independent origins of a unisexual species because of geographic variation in the ancestors. Clones such as C2 and C4 appear to be restricted to a single population, and C3 dominates in one population. Unique clones also dominate in some populations of parthenogenetic D. dahli [37] and D. armeniaca [39]. All populations have C1 and the other clones differ from it by one or two microsatellite repeats only. Thus, the presence of these clones is more likely due to postformation mutations [46], limited dispersal, and genetic drift. Owing to its widespread and ubiquitous distribution, C1 is likely ancestral in D. unisexualis (Table 1) [24]. All other clones have restricted geographic distributions.
Identification of the original area of hybridization can lead to insights and assessments of dispersal, especially  [47,48]. Notwithstanding, D. raddei nairensis occurs sympatrically with D. unisexualis on the western margin of Lake Sevan [30]. This suggests two scenarios for the origin of D. unisexualis. First, the initial C1 arose in the Kuchak region, and then these lizards dispersed eastwardly to other regions (Artavaz, Lchap, Noratus, and other populations). Alternatively, the population at Sevan may have dispersed to western and southern areas. The site of origin remains uncertain, especially since it was proposed to have occurred on the slopes of Mount Aragats [34]. Parthenospecies of Darevskia appear to have evolved recently [10]. Relative to the parental species, they exhibit great mtDNA similarity and low levels of intraspecific variation. Darevskia unisexualis may have originated about 5000 years ago [9], or along with other parthenospecies approximately 200,000-70,000 years ago [34]. Regardless, dispersal resulted in widespread distributions involving many ecological niches.
Because it is not possible to root the network with an outgroup, statistical parsimony network (Fig. 2) has no evolutionary direction. Accordingly, we cannot be certain about the identity the primitive allele. Further, the implied reticulation results in most alleles having equal likelihoods of association with others. These are not inconsequential concerns [49]. The only seemingly unquestionable association is C8 being derived from C2; most other associations remain possible. Regardless, the most widespread allele is consistent with C1 being ancestral.
In summary, microsatellite genotyping analyses [37-39, this study] suggest that clonal diversity in parthenogenetic D. unisexualis and D. rostombekowi, which originated via a single hybridization event, owes to mutations in the initial clones. Similarly, post-formation mutations add to diversity in D. dahli and D. armeniaca, both of which originated via a few hybridizations.

Conclusion
Analyses of four microsatellite loci and single nucleotide polymorphisms (SNPs) in their flanking regions reveal 12 presumptive clones in parthenogenetic D. unisexualis, including one widespread common and 11 rare clones. Assessments confirm that formation of the parthenospecies resulted from the hybridization of female D. raddei nairensis and male D. valentini. Several overall rare clones are numerous and dominate in some populations. Clonal diversity in D. unisexualis appears to result from microsatellite mutations in the initial clone. Parentspecific microsatellite and SNP markers identify multiple clones that allozymes could not. This approach should prove to be equally applicable to detailing the origin and variation of other unisexual species.

Methods
DNA samples were taken from seven populations of parthenogenetic D. unisexualis (n = 109). Analyses also included its parental species: six populations of matrilineal D. raddei nairensis (n = 45), and four populations of D. valentini (n = 17). All samples were from Armenia (Table 3 and Fig. 3).
We used the tips of tails of museum specimens in the herpetological collection of Yerevan State University, as well as a few blood samples collected in 2018 (Table 3). Yerevan State University approved all work with the lizards, which adhered strictly to ethical guidelines. Blood samples were obtained by removing tail tips, which autotomize, and the lizards were then released at the site of collection. DNA extraction used the standard phenol-chloroform method with proteinase K, and resuspension in TE buffer, pH 8.0.
A GenePak PCR Core Kit (Isogene) was used for amplifications in a 20 μl reaction volume, which included approximately 50 ng of DNA and 1 μM of each primer. PCR amplification conditions were used as described previously [39]. Both allelic PCR products of a locus were visualized by electrophoresis in 8% native (nondenaturating) polyacrylamide gel and then excised, purified and sequenced in both directions as previously described [39].
As before [54], a statistical parsimony haplotype network was calculated using TCS v.1.21 to visualize geographic distribution of clones and overall similarity. Notwithstanding, homologous alleles in parthenogenetic clones had linear arrangements via repeat number with little or no recombination. Thus, our coding (Table 4) considered gaps as a fifth state [37,55].
Additional file 2: Table S2. Population indices of gene diversity for four loci in four sampled populations of D. raddei nairensis.
Additional file 3: Table S3. Ethics approval and consent to participate Our study relied mostly on previously collected samples. Thus, permission to collect the DNA samples used in this study was not required. We used tailtips of lizards preserved as museum specimens at Yerevan State University, which were collected in 2001-2005 with capture permit Code 5/22.1/51043 issued by the Ministry of Nature Protection of the Republic of Armenia for scientific studies. We also used blood samples of lizards, collected in 2018 and then these lizards were released in the site of collection, which precluded needing a collecting permit. Sampling was approved by the ethical guidelines of Yerevan State University and did not harm the lizards. The study was carried out in strict accordance with the guidelines and scientific standards. 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.