Parthenogenesis in Darevskia lizards: A rare outcome of common hybridization, not a common outcome of rare hybridization.

Hybridization is a common evolutionary process with multiple possible outcomes. In vertebrates, interspecific hybridization has repeatedly generated parthenogenetic hybrid species. However, it is unknown whether the generation of parthenogenetic hybrids is a rare outcome of frequent hybridization between sexual species within a genus or the typical outcome of rare hybridization events. Darevskia is a genus of rock lizards with both hybrid parthenogenetic and sexual species. Using capture sequencing, we estimate phylogenetic relationships and gene flow among the sexual species, to determine how introgressive hybridization relates to the origins of parthenogenetic hybrids. We find evidence for widespread hybridization with gene flow, both between recently diverged species and deep branches. Surprisingly, we find no signal of gene flow between parental species of the parthenogenetic hybrids, suggesting that the parental pairs were either reproductively or geographically isolated early in their divergence. The generation of parthenogenetic hybrids in Darevskia is, then, a rare outcome of the total occurrence of hybridization within the genus, but the typical outcome when specific species pairs hybridize. Our results question the conventional view that parthenogenetic lineages are generated by hybridization in a window of divergence. Instead, they suggest that some lineages possess specific properties that underpin successful parthenogenetic reproduction.

The role attributed to hybridization in evolution has changed markedly in recent decades. When the biological species concept was first developed (Dobzhansky 1940; Mayr 1942), hybridization between different populations was considered mostly ineffective. Hybrid zones were (and still are, in some cases) expected to be transient and interpreted as evidence of reproductive barriers between the diverging groups (Soltis 2013; Abbott et al. 2016). However, with the development of high‐throughput sequencing techniques, extensive evidence of ongoing hybridization within and between diverging clades has accumulated (Abbott et al. 2016). These techniques enabled the investigation of the role of hybridization in eroding barriers to gene exchange during speciation (Abbott et al. 2013), in accelerating speciation by adaptive introgression (Pereira et al. 2014; Stankowski and Streisfeld 2015; Cronk and Yang 2016), and in initiating the formation of new species (Mavárez et al. 2006; Taylor and Larson 2019).

That interspecific hybridization can play a part in the origin of new sexual species has now been extensively documented (e.g., Mavárez et al. 2006; Pereira et al. 2014; Wang et al. 2017; Schumer et al. 2018; Cuevas et al. 2021). However, an alternative consequence of interspecific hybridization is the origin of asexual forms (Schultz 1973; Vrijenhoek 1989). These forms are able to escape the constraints of sex and, if parthenogenetic reproduction is successful, develop into new asexual “species” (sensu Barraclough 2010). In vertebrates, interspecific hybridization is the process responsible for the origin of all known parthenogenetic species (e.g., Murphy et al. 2000; Alves et al. 2001; Fujita and Moritz 2009; Freitas et al. 2016a; Brunes et al. 2019; Fujita et al. 2020). However, it is not known how the hybrid origin of parthenogenesis relates to the overall history of introgressive hybridization within a given clade. For example, does hybridization between two specific sexual species always generate parthenogenetic descendants, or can it result in introgression as well? Also, it is still unclear whether the origin of vertebrate parthenogenesis is induced by hybridization between species at a specific degree of divergence alone. In other words, whether all species pairs with a specific degree of divergence could generate parthenogenetic hybrid descendants (Balance Hypothesis: Wetherington et al. 1987; Moritz et al. 1989), or only specific phylogenetic clades can act as parental species (Phylogenetic Constraint Hypothesis: Darevsky 1967; Kupriyanova 1997, 1999; Avise 2008).

According to the Balance Hypothesis, parthenogenetic vertebrates arise by the hybridization of two sexual species divergent enough to disrupt meiosis in the hybrids, yet not so divergent as to seriously compromise hybrid viability or (parthenogenetic) fertility (Kearney et al. 2009). Under the Phylogenetic Constraint Hypothesis, asexual lineages arise only from the hybridization between specific sexual species (or clades) with genetic factors that allow them to interbreed and produce viable hybrids capable of reproducing parthenogenetically. When hybridization leads to the formation of fertile sexual hybrids, these hybrids can potentially vector gene flow between the parental species. In contrast, if when two species meet they can only produce parthenogenetic hybrids, then there is no possibility for gene flow among the parental pairs. Therefore, patterns of gene flow between the parental lineages over time can be used to clarify whether distance alone can explain the origin of parthenogenetic hybrids, or if other factors need to be taken into account.

Here we investigated the rock lizard genus Darevskia, a fascinating group with over 25 sexual and seven parthenogenetic species currently described (Arribas 1999). All seven parthenogenetic species have originated by interspecific hybridization between sexual species during the Pleistocene ice ages (Freitas et al. 2016a), but only four sexual species participated in these hybridization events (Darevsky 1967; Murphy et al. 2000). Of these, two always acted as the mother (D. raddei [Boettger 1892] and D. mixta‐1 [Mehely 1909]; Freitas et al. 2019), whereas the other two acted as the father (D. portschinskii [Kessler 1878] and D. valentini [Boettger 1892]), each pair producing more than one parthenogenetic hybrid. According to the currently available phylogeny (Murphy et al. 2000; Murtskhvaladze et al. 2020), the genus contains three main clades (Saxicola, Caucasica, and Rudis), the maternal species belonging to one clade (Caucasica) and the paternal species to another (Rudis). Besides the events that created the parthenogenetic hybrids, there is evidence for ongoing hybridization: between sexual species of the clade Rudis (Tarkhnishvili et al. 2017), backcrosses between parthenogenetic species and their parental species that result in polyploids with decreased fertility (Danielyan et al. 2008; Freitas et al. 2019), and more broadly between sexual species with overlapping distributions (Darevsky 1967).

Specifically, in this study we aimed to address the following questions: (1) Is there gene flow among Darevskia sexual species, and if so, how widespread are the gene flow events? (2) If there is restriction of parthenogenesis to particular clades (or species), is this because hybridization is most common in those clades, or does hybridization occur throughout the genus with only specific species pairs being able to generate parthenogenetic hybrids? (3) Is the result of the hybridization between a given pair of sexual species only parthenogenetic hybrids, or did it also lead to gene flow between them? To answer these questions, we used targeted resequencing of transcriptome‐derived markers to analyze over 600 loci in a phylogenetic framework. A species tree was built to determine the topology and estimate divergence times between species. Introgressive hybridization events were detected using discordance in tree topology, and formal tests for gene flow. We find widespread introgressive hybridization events, both in deep clades and between more recently diverged species. Surprisingly, we find no evidence for gene flow among the parental pairs of the parthenogenetic hybrids. As the gene flow pattern shows, there was widespread hybridization in the genus in the past, so we propose that the origin of the parthenogenetic hybrids is a rare outcome of the total occurrence of hybridization within the genus, but the typical outcome when specific species pairs hybridize.

Material and Methods

SPECIES SELECTION

We included 69 individuals belonging to 14 sexual species in this study (Fig. 1; Table S1). The six parthenogenetic species were not included in our main analyses because our aim was to quantify rates of gene flow between sexual species. However, we confirmed previous inferences of parentage of the six parthenogenetic species using data from 52 parthenogenetic individuals (File S1; Fig. S1; Table S2). All individuals were captured in the field, where they were photographed and had basic measurements recorded for species identification. The tail tip was removed from each individual (by applying pressure with our fingertips), and preserved in ethanol for DNA extraction. Two individuals of the closely related Iranolacerta brandtii (De Filippi,1863) were used as an outgroup (Mendes et al. 2016; Garcia‐Porta et al. 2019). Total genomic DNA was extracted from approximately 30 mg of the tail tip following standard high‐salt protocols (Sambrook and Russell 2001) or the Qiagen DNeasy Blood and Tissue kit.

Figure 1

Range distributions of the sexual Darevskia species analyzed here, after IUCN (2020). Given the high level of overlap between ranges, species were divided into different panels according to phylogenetic clade (see Results). Points correspond to the sample locations (GPS coordinates in Table S1).

CAPTURE SEQUENCE AND GENOTYPING

Probes for the sequence capture were designed from Darevskia transcriptomes (see Supporting Information for details). RNAseq libraries for the four samples were prepared in Edinburgh Genomics, Edinburgh, UK where they were then sequenced in a single lane, using an Illumina HiSeq 2000 machine (100 bp paired‐end reads; insert size around 80 bp). The transcriptomes of two Darevskia species (D. valentini and D. unisexualis [Darevsky 1966]) were de novo assembled into contigs, which corresponded to putative genes. After eliminating contigs duplicated between transcriptomes, they were merged and treated as a single transcriptome reference. The contigs were then filtered and used as a reference for mapping the same RNA reads, subsequently calling SNPs and assessing genetic diversity within the assembled contigs. The most variable contigs (hereafter locus/loci), that is, the ones with most SNPs, were then selected for probe design. Probes were designed so that they were 120 bp in length, and tiled for a 3× coverage, so that each subsequent probe included the last 80 bp of the previous probe. A total of 10,000 probes were developed, spanning 625 loci (corresponding to 575,561 bp in total). Library preparation and Illumina paired‐end sequencing of the captured fragments was performed commercially at Rapid Genomics Inc (Florida 32601 USA), resulting in 100 bp reads. Given the sequencing method used and the good DNA quality of the samples, allelic dropout tests were not considered necessary (see Bose et al. 2018).

Capture sequenced reads were filtered with Cutadapt (Martin 2011) eliminating all reads shorter than 40 bp and with PHRED quality lower than 20. Stampy (Lunter and Goodson 2011) was used to map the reads against the reference contigs used to design the probes. In the case of multiple mappings that are almost (but not exactly) the same, Stampy uses the substitution rate (SR) parameter to choose the best mapping by influencing the mapping quality (MAPQ) score. This parameter was set to take into account the distance of the different species to the reference transcriptome: higher for phylogenetically distant species and lower for phylogenetically close species. Here, D. valentini (a member of the Rudis group) was used as the reference species because the contigs derived from the transcriptome of the hybrid species D. unisexualis were a consensus between D. raddei and D. valentini haplotypes (its parental species). Individuals were divided into three groups: Rudis group samples were mapped with SR = 0.01, non‐Rudis group samples were mapped with SR = 0.03, and the outgroup with SR = 0.1 (note that the mapping was repeated after correcting erroneous morphological identification of some individuals, see below). SR values were selected after trial‐and‐error tests for highest mapping efficiency when comparing numbers of variants (after SNP calling) and flagstat values for the alignments (G. Lunter, pers. commun.). Samtools (Li et al. 2009) was used to convert to bam, sort, index, and calculate basic statistics of the alignment (flagstat option). The MarkDuplicates option in Picard tools was used to remove PCR duplicates, which can be common in a capture sequence protocol. SNP calling was performed with FreeBayes (Garrison and Marth 2012) using minimum coverage of two per position. To account for the expected high divergence between samples (given our analysis is at a genus level, and not species level), we called the four best alleles (instead of the default 2). Post‐vcf filtering was performed with VCFtools (Danecek et al. 2011) and vcflib (https://github.com/vcflib/vcflib), keeping variants with genotype quality >30 and coverage between 5 and 400 for ingroup species (Fig. S2). Ingroup individuals D. parvula (Lantz and Cyren 1913) (sample no: 12738), D. valentini (12641), and D. valentini (16280) were filtered separately because they had higher mean coverage; the selected coverage interval was 5–800 read depth. Outgroup individuals (I. brandtii) were filtered with the same parameters as the ingroup except for coverage (read depth = 2–400). Only SNPs were used in the downstream analyses and indels were removed. Multiple nucleotide positions (MNP) identified by FreeBayes were deconstructed with the vcfallelicprimitives option from vcflib.

SPECIES STRUCTURE AND POPULATION ANALYSES

To identify any individuals that were misidentified on the basis of morphology, a discriminant analysis of principal components (DAPC) was performed on the SNP dataset (biallelic loci only) using the ADEGENET package (Jombart et al. 2010). The optimum number of clusters was chosen with the find.clusters option, by comparing the different clustering solutions using a Bayesian Information Criterion (BIC). Genetic composition was further investigated with STRUCTURE (Pritchard et al. 2000) on the total dataset (not only bi‐allelic SNPs), which also allowed us to estimate the level of admixture for each species. A random SNP per locus was selected and 10 runs were performed, each with a different combination of randomly selected SNPs. Values from the different runs were combined with Clumpp (Jakobsson and Rosenberg 2007). To select for the best number of clusters K, we used Delta K (Evanno et al. 2005) and the ADEGENET best K approximation calculated with the BIC. The K value with the highest concordance between methods was chosen.

PHYLOGENETIC ANALYSES

Alignments of all samples for each marker were made from the vcf file using a custom python script (https://github.com/SusanaFreitas/capture‐data). These alignments were then trimmed to remove missing positions allowing a maximum of 10% missing genotypes for each position in each alignment, because missing data can have a strong negative impact on phylogenetic analyses. Each individual marker was analyzed with Bayesian inference (BI) (MrBayes version 3.2; Huelsenbeck and Ronquist 2001; Ronquist and Huelsenbeck 2003) and the resultant output trees were later used as input for the Bucky analysis (Larget et al. 2010; see below for detailed information). In MrBayes, analyses started with randomly generated trees and ran for 30 × 106 generations, with sampling at intervals of 1000 generations, producing 30,000 trees. A GTR + I + Γ model (a General Time Reversible model with invariant sites and a gamma‐shaped distribution of rates across sites) was chosen because the data were sufficient to estimate all parameters. Two independent runs were performed on each dataset. To ensure the analysis was performed on the stabilized likelihood, 40% of each total run was treated as burn‐in.

Phylogenetic relationships were constructed with *BEAST2 version 2.4.4 (Bouckaert et al. 2014), and the species tree topology was then compared with topologies calculated by SVDQuartets (Chifman and Kubatko 2014; Kubatko and Chifman 2019) and the primary concordance tree from Bucky, the latter also being a reliable estimate of the species tree under the coalescent model (Degnan and Rosenberg 2009) (for more info, see below). The *BEAST2 package performs simultaneous gene tree and species tree inferences under a combined coalescent and sequence evolution model. To optimize the time and resource consumption in the *BEAST2 analysis, we selected a marker subset consisting of the 31 most informative markers (which corresponded to a combination between the longest markers and the ones with the highest number of SNPs) (Ogilvie et al. 2016). Tests of positive and purifying selection (codon‐based Z‐test) were performed on two markers (MC1R and C‐mos), for which we had confidently assigned reading frames, using MEGA (Kumar et al. 2018) and it showed no evidence for deviation from neutrality. The tMRCA was estimated setting as a prior the split between the ingroup (Darevskia) and the outgroup (Iranolacerta) previously estimated as 10 million years ago (Mendes et al. 2016), and run under a strict clock prior. To improve the convergence of the runs, this split was calibrated ignoring uncertainty. To ensure run convergence for *BEAST2, we required ESS values above 200 for all model parameters. Posterior tree samples were combined to produce maximum clade credibility species trees with Tree Annotator version 2.4.4 (Bouckaert et al. 2014), with node heights set to mean age estimates. SVDQuartets (Chifman and Kubatko 2014, 2015) was used with the most variable 300 markers concatenated and run together. This method calculates the unrooted topology for quartets of taxa for each value of a matrix of site patterns (the matrix of concatenated SNPs), and later infers the species phylogeny using the quartet methods (Chifman and Kubatko 2015). SVDQuartets does not allow for clock rate assumptions, so the resulting tree branch lengths were not estimated.

TESTS OF INTROGRESSION

Bucky (Larget et al. 2010) uses Bayesian concordance analysis to estimate the general phylogenetic history of a group of taxa while taking into account the incongruence among gene trees. This software does not make any assumptions regarding the reasons for discordance, which could be either incomplete lineage sorting (ILS), introgressive hybridization, or horizontal gene transfer (HGT); downstream comparative analyses should be undertaken to understand and distinguish the causes of discordance. Bucky takes as input Bayesian inference (BI) trees estimated for each locus analyzed. The higher the number of taxa analyzed, the larger the tree space, which reduces the probability of detecting the same tree more than once. Therefore, to increase Bucky's sensitivity, taxon number was reduced by merging some of the most recently diverged species (identified in the *BEAST2 analysis) into 13 terminal clades. The species merged were D. brauneri [Mehely 1909] and D. saxicola [Eversmann 1834]; D. mixta‐2 and D. parvula; and the two D. praticola [Eversmann 1834] clades. Only one randomly selected sequence per taxon was used in the BI runs. BI trees were estimated with MrBayes for each marker, retaining one instance of each topology found. Given the low variability and thus low phylogenetic information of some markers, only the 300 most variable out of the 625 targeted loci were selected for the Bucky run. Next, the total tree distribution for all markers used was examined and, with a nonparametric Dirichlet process, trees were clustered into groups that share the same topology and later used to estimate the concordance factor (CF). The Dirichlet process uses an a priori level of discordance α (here α = 1), which combines the information in the sequences of those genes that are inferred to be congruent (Chung and Ane 2011). CFs for each clade were estimated based on the proportion of trees that contain that clade, and a primary concordance tree was then constructed to include the clades with the highest estimated CFs. Clades not present in the concordance tree but with significantly high CF values were considered discordant.

Patterson's D‐statistics (Durand et al. 2011) was used to test whether the pattern of shared variation between species provided by Bucky was better explained by gene flow than by incomplete lineage sorting. These tests, also called ABBA/BABA tests, determine whether the number of derived alleles shared by reciprocally monophyletic taxa is greater than expected under the null model of incomplete lineage sorting (ILS). Ancestral (A) and derived (B) alleles are identified by comparison with an outgroup taxon. The numbers of ABBA and BABA patterns are calculated in a group of four clades: two sister clades (P1 and P2), a putative admixed clade (P3), and the outgroup (O). If no differential gene flow has happened between one of the sister taxa (P1, P2) and the “admixed” taxon (P3), the numbers of ABBA and BABA patterns are expected to be the same. D‐statistics were estimated with the doAbbababa option in ANGSD (version 0.930) (Korneliussen et al. 2014), which takes as input the bam files of each taxon (P1, P2, P3, and O), and outputs the D‐statistic and a Z‐score for each test. Because D‐statistics have to be estimated respecting the topology of the taxa tested, the number of tests we could perform within the genus was limited. Respecting the phylogenetic relationships observed from our previous analyses (*BEAST2 and Bucky concordance tree), a large majority of all possible taxa combinations were tested for introgression. For each tested combination of taxa (hereafter each “test”), individuals of each species were randomly chosen repeatedly so that, for each test, D‐statistics were calculated for multiple combinations of individuals. Given the lower coverage of the outgroup Iranolacerta, to increase the number of positions tested we used D. valentini as an outgroup whenever possible, which should increase the power to detect gene flow. For each combination of taxa, 150 random combinations of individuals were run. Significance of tests was interpreted as follows: strong evidence for gene flow when mean Z‐score per test was ≥4, good evidence for gene flow when mean Z‐score was between 3 and 4, and suggestive evidence for gene flow when mean Z‐score was <3, but Z was ≥4 for some combinations of individuals. In the standard four‐individual (or taxon) tests, a Z‐score ≥3 (P = 0.01) or Z ≥ 4 (P = 0.001) indicates introgression between P1 and P3 (if the D‐statistic is <0), or between P2 and P3 (if the D‐statistic is >0), compared with the expectation from ILS alone (Korneliussen et al. 2014).

In addition, formal tests for introgression with three‐population (f 3) and four‐population (f 4) statistics (Reich et al. 2009) were used to estimate admixture proportions in groups of populations (or taxa) because these tests differ in sensitivity to admixture under different conditions. f 3 and f 4 statistics estimate correlations between allele frequency differences of three (f 3) or four (f 4) populations joined by an unrooted tree. If the populations tested are not admixed, genetic drift (estimated from the allele frequency as a proxy) should generate uncorrelated allele frequencies. If correlation between populations is found, this is indicative of gene flow. The f 3 statistic (X; Y, W) tests for admixture between a test population X and two reference populations, Y and W. In the case of no admixture, its expected value is positive, whereas it is negative in the case of admixture between X and Y, X and W, or both Y and W with X. On the other hand, the f 4 statistic tests the correlation between allele frequency differences in an unrooted tree (A,B; C,D). It is equal to zero if there are no correlations and no admixture occurs across the tree branches; positive if admixture occurs between A and C, between B and D, or both; and negative if admixture occurs between A and D, B and C, or both. The statistical significance for both tests was assessed using a Z‐score (f 3 or f 4 divided by its standard deviation), and a threshold Z‐score of 3 corresponding to P = 0.01 was used. The programs threepop and fourpop in the TreeMix package (version 1.13; Pickrell and Pritchard 2012) were used to calculate f 3 and f 4 statistics, respectively, as well as Z‐scores to assess the significance of the test values for all possible combinations of three and four taxa within Darevskia. These Z‐scores were then transformed into P‐values and significance was confirmed after an False Discovery Rate correction (Benjamini and Hochberg 1995) for multiple comparisons.

To visualize the extent, directionality and time of gene flow, we used the graph‐based method developed in TreeMix. This approach uses allele frequencies at genome‐wide polymorphisms and a Gaussian approximation of the genetic drift among populations to estimate population graphs where the tree‐like relationships among populations are supplemented by a number of admixture events within a chosen interval, here 0–21. The best number of admixture events was estimated by comparing the proportion of the variance in population covariances explained by each population graph with a different number of admixture events, and the best model was also expected to show a consistent improvement in the percentage of variance explained across replicate runs. Following this, migration edge values were added to the tree, and the estimated admixture events were observed. For TreeMix, f 3 and f 4, a total of 5823 SNPs was used, missing data excluded. Only ingroup species were tested.

Results

MAPPING AND SNP CALLING

We developed and sequenced 625 new markers for our 69 samples from sexual species (Table S1) and 52 parthenogenetic hybrids (Table S2). In our main dataset (with the sexual species only), our sequences comprised 574,341 bases in the transcriptome reference alignment used, corresponding to 429,078 sites recovered with our probes in the 69 individuals analyzed. On average, 226 bp were missing from each locus (or reference contig), and the amount of recovered sequence was between 25 and 2893 bp per target (reference loci ranged from 550 to 3740bp). All loci were variable, ranging from 1 to 379 SNPs in total (0.12%–24% variable positions excluding the outgroup). In total, the dataset had 64,150 SNPs, and of these 45,959 positions showed a minor allele count of 2 or higher excluding the outgroup, that is, they were phylogenetically informative. Diversity (nucleotide diversity, pi, and Tajima's D) among species did not show any systematic differences between species (Figs. S3, S4).

Analyses including the parthenogenetic taxa are reported in File S1. These analyses confirmed the previously described parentage of parthenogenetic species (Fig. S1).

PHYLOGENETIC INFERENCE

Overall, the topology of all the species trees estimated here (for sexual species) recovered the three main clades: Saxicola, Caucasica, and Rudis. The different analyses recovered the same topology, except for the branching of D. parvula. Darevskia parvula and D. mixta‐2 were found to be the same species, and D. mixta‐2 a morphological misidentification of a new population not described before. The D. parvula taxon branched within the Caucasica clade, as a sister species of D. clarkorum/D. mixta‐1, in the *BEAST2 species tree (Fig. 2), but as an outgroup to both the Caucasica and Saxicola clades in Bucky and SVDQuartets trees (Figs. 3, S5). The high level of gene flow found between D. parvula and the Rudis clade (see INTROGRESSION TESTS) likely explains this disagreement. Hence, we used the *BEAST2 topology (positioning D. parvula within the Caucasica clade) (Fig. 2) in the D‐statistics tests. Finally, the topology recovered in our *BEAST2 species tree overall coincided with previous phylogenetic studies on the genus (Murphy et al. 1996; Freitas et al. 2019, 2016b; Murtskhvaladze et al. 2020), but provided greatly increased support for the phylogenetic relationships because of the larger number of markers included.

Figure 2

DensiTree representation of the *BEAST2 species tree inference for the sexual Darevskia species analyzed in this study based on the 31 most informative markers. Parental species of parthenogenetic lineages are highlighted: green for the paternal species and gray for the maternal species. The major clades recovered are noted with letters, (a) Rudis, (b) Caucasica, and (c) Saxicola. Branch colors correspond to the likelihood of the trees: blue is the most likely tree topology, red the second, and green the third. The consensus tree is drawn as a single blue line above all the others and its posterior probability values are presented on each node, and noted as a black dot whenever they are 1.0.

Figure 3

Primary concordance tree inferred with Bucky, constructed with the 300 most variable out of the 625 targeted loci. Insets correspond to frequent clades that are not represented in the primary concordance tree. Vertical lines connecting branches indicate significant incongruences, with colors indicating concordance factor (CF) values as identified in the gradient (gradient bar on the bottom right). Node values correspond to the concordance factor for the respective clade. Branch length measured as coalescent units (scale bar).

Our divergence time estimates suggest that the Darevskia genus started differentiating during the early Pliocene, ∼5 million years ago, splitting into three major clades, Saxicola, Caucasica, and Rudis (Fig. 2). The first node in the estimated phylogenetic tree for Darevskia (Figs. 2, S6) separates the Rudis clade from the remaining taxa, and is dated ∼5 [3.9, 6.1] million years ago (posterior mode [95% credible interval]). The next node separates the Saxicola clade (D. praticola, D. brauneri, and D. saxicola) and the Caucasica clade (D. raddei, D. mixta, D. parvula, and D. chlorogaster [Boulenger 1908]; D. steineri [Eiselt 1995]; D. defilippii [Camerano 1877]; D. clarkorum [Darevsky and Vedmederja 1977]; D. derjugini [Nikolsky 1898]), at ∼4.5 [3.5, 5.5] million years ago. The Balkan and Caucasian populations of D. praticola were estimated to have separated ∼1 [0.651, 1.463] million years ago. However, these time estimates must be interpreted with caution, given the evidence for introgression described below.

INTROGRESSION TESTS

We found evidence of widespread introgression among species within the genus Darevskia, both recent and past. Surprisingly, no gene flow was detected among the parental pairs of the asexual hybrids.

Bucky found 13 incongruences (Fig. 3), eight of which coincided with either significant D‐ or f tests (see Table 1 for details). Most incongruences that did not coincide with the other gene flow tests were incongruences between deep clades and not terminal taxa, making them difficult to detect with the other tests. The only recent Bucky incongruences that did not coincide with either D or f tests, or TreeMix, were the ones between D. raddei and D. steineri, and between D. rudis (von Bedriaga 1886) and D. portschinskii that are likely due to incomplete lineage sorting and not gene flow. In general, the incongruences (all 13, but also the subset coincident with the gene flow tests) were distributed throughout the depth of the tree, from deep branches (e.g., Rudis clade‐D. parvula) to recently diverged taxa (e.g., D. defilippii‐D. raddei).

Table 1

Summary table of gene flow evidence detected

Summary table of all significant instances of evidence for gene flow between terminal branches found by the methods used here. Events involving deeper branches (such as some of those detected by Bucky and TreeMix) are not represented here. Entries in the lower triangle show the significant instances of gene flow detected by each method used: Bucky (1), ABBA/BABA (2), TreeMix (3), and f4 (4). Each significant f4 test is identified by a letter from a to p (4a–4p). The top half of the table indicates which ABBA/BABA tests were significant (see Table S3 for test details). Green squares indicate species pairs that have given rise to parthenogenetic hybrid lineages. Species names are as follows: D. brauneri (BRAU), D. chlorogaster (CHLO), D. defilippii (DEFI), D. derjugini (DER), D. mixta‐1 (MIX1), D. parvula (MIX2 and PAR), D. portschinskii (POR), D. praticola (PRA), D. raddei (RAD), D. rudis (RUD), D. saxicola (SAX), D. steineri (STEI), and D. valentini (VAL).

D‐statistics provided evidence for gene flow in 20 tests, distributed throughout the genus, from a total of 147 tests performed (Table S3). Three tests showed strong evidence for gene flow; between D. chlorogaster and D. praticola, D. clarkorum and D. praticola, and D. mixta‐1 and D. praticola. Eleven tests showed good evidence for gene flow and six other tests showed suggestive evidence (for a summary of test results, see Table 1).

The f4 analysis also provided additional evidence for widespread admixture in the genus. Given the nature of the f4 test, which analyses two species pairs and should be evaluated as admixture between either/or both species pairs tested, it was only used here as a confirmation of other admixture evidence. A total of 16 tests out of all possible combinations analyzed (2575) showed evidence of admixture (see Table 1 for overlap with other evidence for gene flow). The significant f4 test 4h (Table 1) included the pairs D. portschinskii/D. mixta‐1 (parents of the parthenogenetic D. rostombekowi) and D. valentini/D. rudis. Given the nature of the f4 test, its significance means that there is evidence for gene flow either within both pairs or only one of them. Because there was strong evidence for gene flow for the pair D. valentini/D. rudis from other tests (Bucky and five other f4 tests), and no other introgression test detected gene flow between D. portschinskii and D. mixta‐1, it is likely the gene flow between D. valentini and D. rudis was solely responsible for the significance here. No significant values were obtained from the three‐population formal test (f3).

TreeMix was used to construct a tree connecting the 14 species to describe and position the admixture events. The variance explained by the TreeMix ML model without migration edges was 98.9% (Fig. S7), and it was significantly increased by adding migration edges up to the model with four migration edges (Figs. 4, S7). No evidence for introgression between the parthenogens' parental pairs was found by TreeMix (Fig. 4). The TreeMix model also showed admixture events throughout the evolutionary time‐scale of Darevskia, with events happening recently (such as the D. raddei‐D. chlorogaster introgression) and others estimated to predate the first split within its clades (such as the D. praticola‐Caucasica clade). This analysis also showed, in concordance with all the others, introgressive events between the parthenogens’ sexual parentals and other Darevskia species, but not within the parthenogens’ parental pairs. The values of the residual covariance between each species pair identify deviations from the proposed model. For those species pairs with darker colors in the matrix of Figure 4, the bifurcating model together with the estimated migration events may still not fully explain the covariance among Darevskia species.

Figure 4

TreeMix Maximum Likelihood (ML) tree estimated from the allelic frequencies of 5823 SNPs. ML tree (top) and residuals from the fit of the ML model (bottom) inferred with 4 migration edges. The arrows on the ML tree indicate the directionality of gene flow on the migration edge and the color of the edge reflects the intensity of admixture. The heat colors in the residual matrix represent the covariation unexplained by the corresponding tree (light colors are close to 0, dark colors have higher residuals). Positive residuals indicate taxa are more closely similar to each other than expected under the proposed model.

Overall, all introgression tests (Bucky, D‐ and f4 tests, and TreeMix) tended to show consistent patterns (Table 1). Although some species were frequently implicated in gene flow events (for instance D. raddei was found in 11 significant events of gene flow with other species), others showed little evidence of gene flow (such as the case of D. saxicola, which was picked only once by f4 with D. portschinskii). However, and despite the parental species being some of the species more frequently implicated in the gene flow events, all introgression analyses undertaken agreed in the absence of evidence for gene flow among the parental pairs of the parthenogenetic hybrids.

Discussion

In this study, we assessed the extent of introgressive hybridization among 14 sexual species of Darevskia rock lizards. In vertebrates, hybridization has repeatedly generated parthenogenetic hybrid lineages. However, it is not known whether this hybridization is a specific property of the parthenogens’ parental species pair, nor if the result of the hybridization was only hybrid asexuality, or perhaps also gene flow between the parental species. In this study, we tested for evidence of gene flow between sexual species to understand the relationship between hybrid asexuality and reticulate evolution. We found widespread gene flow, happening throughout the divergence continuum, both between deep branches and following more recent splits, but not between the parental pairs of the parthenogenetic hybrids, showing that the origin of parthenogenesis in Darevskia is a rare outcome of common hybridization within the genus, but also a common outcome of rare hybridization between specific pairs of species.

EVOLUTIONARY HISTORY OF Darevskia

Our species tree is consistent with earlier Darevskia phylogenetic studies (Murphy et al. 2000; Murtskhvaladze et al. 2020), except for the placement of two species (D. parvula and D. mixta‐1, hereafter D. mixta). Earlier phylogenies were based on a limited set of markers (e.g., Murphy et al. 2000; Ahmadzadeh et al. 2013; Pyron et al. 2013; Roquet et al. 2014; Murtskhvaladze et al. 2020) and/or fewer species, excluding several of the Darevskia species found only in Iran (currently, D. chlorogaster and D. defilippii species complexes with D. steineri as part of the latter; Ahmadzadeh et al. 2013), and the Balkan clade of D. praticola, all of which we include here. In our study, D. parvula is shown to be an early diverging lineage within the Caucasica clade, whereas other studies placed D. parvula within the Rudis clade (Fu et al. 1997; Murphy et al. 2000), or as an outgroup to the other Darevskia species (Garcia‐Porta et al. 2019; Murtskhvaladze et al. 2020). This positioning is likely a misplacement due to the extensive gene flow we found in this study between D. parvula and the Rudis clade (Figs. 3, 4; Table 1). Similarly, in our analyses D. mixta is consistently found together with D. clarkorum, but reported as a sister species to D. raddei in previous studies (Murphy et al. 1996). Even though later studies have also recovered D. clarkorum as the sister species to D. mixta (Garcia‐Porta et al. 2019; Murtskhvaladze et al. 2020), the relationships within this clade (the Caucasica clade) are somewhat different to the topology we recover. Like D. parvula, D. mixta, D. raddei, and D. chlorogaster (all from the Caucasica clade) also showed evidence for frequent introgressive hybridization with other sexual species (Fig. 4; Table 1), which could explain why it is difficult to place them correctly without using an approach that takes into account the influence that strong introgression can have on species tree estimation. Most previous Darevskia phylogenies are based on concatenation methods (Pyron et al. 2013; Roquet et al. 2014; Garcia‐Porta et al. 2019; Murtskhvaladze et al. 2020), which cannot handle heterogeneity among gene trees and have been consistently shown to perform poorly when gene flow (or other evolutionary processes) causes incongruence in the gene trees (e.g., Liu et al. 2015; Jiang et al. 2020). On the other hand, the multispecies coalescent methods used here have been shown to perform well in estimating species trees even in the presence of gene flow (SVDQuartets: Long and Kubatko 2018). These methods can integrate gene tree incongruences (caused by incomplete lineage sorting or gene flow) into the species tree construction (Bucky: Larget et al. 2010), or co‐estimate gene trees with the species tree (Heled and Drummond 2010; Bouckaert et al. 2014; Mallo and Posada 2016; Ogilvie et al. 2016).

The first split within the genus Darevskia happened in the early Pliocene, ∼5 million years ago. This period was characterized by high volcanic activity in the study area (that lasted until the Mid Pliocene: Tarkhnishvili et al. 2013). This volcanic activity contributed to great orogenic and climate changes that led to a regression of seawater and the formation of the modern outlines of the Black and Caspian Seas among other consequences (Popov et al. 2006). These abrupt changes offered new habitats and allowed for diverging taxa to expand into new territories, leading to rapid divergence and paving the way for the high speciation rates of this genus. We also find that the divergence within D. praticola, between the Balkan and Caucasian clades, is not as old as previously reported (2.5 million years: Freitas et al. 2016b), but instead is estimated to have occurred ∼1 million years [0.651, 1.463] ago. This contradiction further confirms the idea that multilocus coalescent phylogenies, which take into account different gene trees, are better suited to get a more accurate evolutionary history of groups of organisms with such complex histories of reticulate evolution and hybridization.

HYBRIDIZATION AND GENE FLOW

Although we cannot formally date the introgressive events on the tree, the different introgression analyses performed here can detect gene flow that happened in different periods. The formal tests (D‐ and f4 statistics) are not sensitive to recombination because they can detect gene flow even after recombination has partially diluted the introgressed genome. In contrast, because Bucky is based on gene tree analyses, it can be more affected by recombination. If recombination breaks down the introgressed sequence analyzed, clades may not coalesce, thus diluting the signal of the introgression. Following this reasoning, although D‐statistics and f4 formal tests can detect both old and recent introgressive hybridization, Bucky can only detect gene flow events that have not been diluted by recombination, and therefore are relatively recent. A caveat of Bucky analyses is that it does not distinguish between ILS and gene flow (Larget et al. 2010). However, because the formal tests (D and f4 statistics) used are not expected to be affected by ILS (when significant, gene flow is the only responsible actor; Durand et al. 2011), they were used to confirm whether incongruences detected by Bucky were due to ILS or gene flow. To complement these analyses, TreeMix was used not only to confirm the patterns of introgression found, but also to further place the historical splits and admixture within Darevskia. A final caveat in relation to the gene flow tests performed here: because we use protein coding sequences (and not neutral markers), we have to consider that these markers may be under purifying selection, which would tend to cause underestimation of introgression. However, the markers tested showed no evidence of being under selection (despite the same markers being found to be under selection in other reptile groups). Moreover, it is difficult to see how possible purifying selection would generate a bias against detecting introgression between parents of asexual taxa. Therefore, we think putative purifying selection is unlikely to have influenced our major results.

All tests of introgression used here showed evidence for gene flow events, except the f3 formal test. A significantly negative f3 (X; Y, W) value indicates that the frequencies of alleles in population X tend to be intermediate between Y and W, which can only arise if population X resulted from a mixture of populations related to Y and W. However, if populations had had enough time for drift, or X had suffered bottlenecks after the admixture event, consistent with the complex demographic history of this genus, the f3 statistic may not recover a signal of introgression (Decker et al. 2014; Kamdem et al. 2017), which can explain why we find evidence for gene flow with the other tests, and not the f3.

We find widespread gene flow events between taxa and along the divergence continuum, not only restricted to those recently separated, but also between those more anciently diverged, even though the precise date and extent of introgression could not be measured. Given the reticulate evolutionary history of this genus, the parthenogenetic hybrid species are present witnesses to past interspecies hybridization events. It has been suggested that secondary contact with hybridization was facilitated by the complex landscape in the Anatolia/Caucasus region. Indeed, secondary contact and introgressive hybridization have happened frequently not only in this genus (Tarkhnishvili et al. 2013; Freitas et al. 2016a; Candan et al. 2021; Rato et al. 2021), but also in many other vertebrate genera occupying the same region (e.g., Pabijan et al. 2017; van Riemsdijk et al. 2017; Wielstra et al. 2017; Levin et al. 2019). The evidence for widespread introgressive hybridization in Darevskia confirms the earlier expectations for frequent secondary contact throughout the evolution of this group.

GENE FLOW, PARTHENOGENESIS AND REPRODUCTIVE BARRIERS IN DIVERGING CLADES

Widespread gene flow, as we show in Darevskia, has been repeatedly found in many other vertebrate groups, including fish (cichlids: Malinsky et al. 2018) and other reptile groups, particularly lizards (Gottscho et al. 2017; Grummer et al. 2018; Caeiro‐Dias et al. 2021). However, most groups with hybridization did not generate parthenogenetic species (at least none that survived to the present). Thus, parthenogenetic vertebrates originate, or persist, in groups with widespread introgressive hybridization (our study; Barley et al. 2019; Morando et al. 2020), but widespread hybridization alone is not sufficient for generating asexual hybrid species.

In Darevskia, the restriction of parthenogenesis to crosses between specific sexual species is not because hybridization is most common between these species, but instead because hybridization is widespread yet only generates parthenogenesis in some clades. Indeed, despite the high number of introgressive events detected among Darevskia in general, no gene flow was detected among the parental pairs of the parthenogenetic hybrids, D. raddei or D. mixta (maternal species) and D. portschinskii or D. rudis/D. valentini (paternal species). This pattern, hybridization events responsible for the origin of parthenogenetic hybrids that cannot produce sexual (and fertile) hybrids simultaneously, has been described previously in taxa that generated parthenogenetic individuals via experimental crosses between two sexual species (fish: Schultz 1973; nematodes: Lamelza et al. 2019). Furthermore, other lizard groups with sexual and parthenogenetic hybrids also present the same pattern of widespread introgressive hybridization, but absence of gene flow between the parental pairs of asexual species (Aspidoscelis: Barley et al. 2019; Liolaemus: Morando et al. 2020).

Considering the evidence for widespread gene flow we found in Darevskia, together with the absence of gene flow among the parental pairs of the parthenogens, we can suggest that these sexual species pairs became isolated early in their divergence, either geographically or by intrinsic barriers, and when hybridizing they could only produce parthenogenetic (and perhaps sexually infertile) hybrids. If the origin of parthenogenetic hybrids in vertebrates was due to the phylogenetic distance between parentals alone, as the Balance Hypothesis suggests, diverging taxa would remain able to exchange genetic material (when in contact) until they reached a certain genetic distance. Then, other groups with widespread introgressive hybridization would be likely candidates to contain parthenogenetic species, because there would have been multiple opportunities to generate them at the appropriate stages of divergence. Moreover, there would be no reason why the parental pairs would not be able to exchange genes (if meeting before reaching the genetic distance that would prevent production of viable and fertile sexual hybrids, but allow them to produce a parthenogenetic hybrid). Given only a subset of species pairs can produce parthenogenetic hybrids, and these species have shown a history of frequent past introgressive hybridization with other sexual species, but not among themselves, the overall evidence is more concordant with the expectations of the Phylogenetic Constraint Hypothesis. This hypothesis suggests that the potential for an interspecies hybridization event to generate a parthenogenetic hybrid is dependent on lineage‐specific genetic factors. Obviously, a certain distance between the parentals is needed to disrupt the sexual viability of the hybrid and hybridization at this early stage might result in gene flow. However, distance alone cannot explain the origin of the asexual hybrids, and intrinsic aspects of the parental pairs must be involved.

Sexual reproduction is widespread in nature, and its higher potential for adaptation is generally considered to explain its prevalence over asexual reproduction (McDonald et al. 2016; Luijckx et al. 2017). However, the demographic advantage of asexual populations (Maynard Smith 1971), the existence of lineages known to be reproducing solely asexually for millions of years (Judson and Normark 1996; but see Schwander 2016), the lack of associated costs to asexuality (in the specific case of Darevskia: Arakelyan et al. 2019), and the potential of asexuals to occupy more diverse habitats than their sexual relatives (in Darevskia: Tarkhnishvili et al. 2010) indicate that a complementary explanation for the scarcity of asexual lineages is needed. If the conditions that lead to the origin of asexual lineages are very particular, this could be a major restriction acting on the source of new asexual individuals and another factor to consider regarding the ubiquity of sexual reproduction in nature (Vrijenhoek 1989; Schwander and Crespi 2009). Our results suggest that parthenogenetic hybrid lineages can only be generated under rather special conditions, via crosses between specific parental species, thus constraining the number of parthenogenetic vertebrate lineages on the basis of a low origin rate.

Our results challenge the conventional understanding of the conditions required for the origin of parthenogenetic lineages by interspecific hybridization, and suggest that genetic distance alone is not a sufficient explanation for hybrid parthenogenesis. Hence, specific properties of the parental genomes and incompatibilities arising in the hybrid genome should be considered in the discussion surrounding the mechanisms behind vertebrate parthenogenesis. Other authors (Kupriyanova 1997, 1999; Murphy et al. 2000) have suggested that the level of divergence must have acted together with sex chromosome incompatibilities in the hybrid. The relevance of sex chromosomes in the mechanism behind asexual reproduction in Darevskia has been further supported by recent evidence that has shown that meiotic cells in D. unisexualis, one of the most widespread parthenogenetic Darevskia, have normal synaptic complexes between all autosomes, but no synapsis is found between sex chromosomes (Spangenberg et al. 2020). Even though lack of synapsis in other vertebrate groups leads to arrest of meiosis and apoptosis of the dividing cells (Bhattacharyya et al. 2013), in rock lizards (or at least in D. unisexualis parthenogenetic hybrids) there is a reduction or even complete absence of these checkpoints, and meiosis proceeds until the formation of oocytes (Spangenberg et al. 2020). Although the lack of checkpoints may be relevant for the mechanisms of parthenogenesis in Darevskia, in other lizard groups, such as Aspidoscelis neomexicana, checkpoints have to be respected for the occurrence of meiosis, and only cells that have had a premeiotic chromosome duplication, with pairing and recombination occurring between identical chromosomes (instead of homologs like in sexual meiosis), are able to finish meiosis and reach the oocyte stage (Newton et al. 2016). Unlike Aspidoscelis, Darevskia parthenogenetic meiosis does not have (or require) a premeiotic chromosome duplication. This disparity in only two groups analyzed suggests that the mechanisms of meiosis in vertebrate parthenogens may not conform to a one‐size‐fits‐all scenario. Indeed, understanding how these organisms escape the constraints of sexual mechanisms for reproduction is important not only to shed light on asexual reproduction itself, but will also help us understanding key mechanisms underlying sexual reproduction and meiosis.

AUTHOR CONTRIBUTIONS

SF, RKB, MAC, and DJH designed the study. SF, MA, ÇI, YK, MAC, and DJH collected the samples. SF generated the molecular data. SF carried out all analyses with input from RKB and AW, and interpretations were discussed between all authors. SF wrote the draft manuscript with all authors contributing to the final version.

Associate Editor: A. G. McAdam

Handling Editor: A. G. McAdam