Allopolyploidy: An Underestimated Driver in Juniperus Evolution

Allopolyploidy is considered as a principal driver that shaped angiosperms’ evolution in terms of diversification and speciation. Despite the unexpected high frequency of polyploidy that was recently discovered in the coniferous genus Juniperus, little is known about the origin of these polyploid taxa. Here, we conducted the first study devoted to deciphering the origin of the only hexaploid taxon in Juniperus along with four of its closely related tetraploid taxa using AFLP markers with four primers combinations. Phylogenetic analysis revealed that the 10 studied species belong to 2 major clusters. J. foetidissima appeared to be more related to J. thurifera, J. sabina, and J. chinensis. The Bayesian clustering analysis showing a slight variation in genetic admixture between the studied populations of J. foetidissima, suggesting an allopolyploid origin of this species involving J. thurifera and J. sabina lineages, although a purely autopolyploidy origin of both J. thurifera and J. foetidissima cannot be ruled out. The admixed genetic pattern revealed for J. seravschanica showed that the tetraploid cytotypes of this species originated from allopolyploidy, whereas no clear evidence of hybridization in the origin of the tetraploid J. thurifera and J. chinensis was detected. This study provides first insights into the polyploidy origin of the Sabina section and highlights the potential implication of allopolyploidy in the evolution of the genus Juniperus. Further analyses are needed for a more in-depth understanding of the evolutionary scenarios that produced the observed genetic patterns.


Introduction
Natural hybridization and polyploidy have long been considered as partners in plant evolution and promoters of biodiversity [1]. Allopolyploidy is when genome duplication involves the genomes from diverged species through hybridization, inducing diversification and speciation. Indeed, allopolyploidy is a common pathway for hybrid speciation [2][3][4]. It has been well documented in angiosperms, yet overlooked in gymnosperms (except Ephedra L. [5,6]). Recently, an unexpected frequency of polyploid species (15% of natural taxa) was determined in the coniferous genus Juniperus L. It was suggested that at least ten polyploidy events have occurred in the course of Juniperus evolution [7]. All polyploid Juniperus taxa are distributed in the old world and mainly in the Mediterranean region [7]. Table 1. Details on the studied populations of Juniperus species. *, material collected and used in the Ph.D. research of Bouchra Douaihy [25]. The ploidy level of each taxon was extracted from Farhat et al. [7]. The ploidy level appended by "i", "ii", or "iii" is the exact samples measured by Farhat et al. [7,12,23].

Amplified Fragment Length Polymorphism (AFLP) Analysis
The AFLP method used is based on the standard protocol [27]. To improve the reliability of AFLP data, 34 replicates of AFLP profiles were generated from independent restriction digests of the same DNA (repetition type 1) and from different DNA extractions of the same sample (repetition type 2) (  (Table 3). Polymerase chain reaction (PCR) was carried out in a final volume of 50 µL. In each reaction, 5 µL of diluted restrictionligation products was added to 10 pmol of each primer, 0.2 mM of each dNTP, 2.5 mM of MgCl 2 , 1X Taq DNA polymerase buffer without MgCl 2 , and 1 unit of Taq DNA polymerase (MP Biomedicals, Santa Ana, CA, USA). The pre-selective amplification protocol involved initial denaturation at 94 • C for 1 min, followed by 20 cycles at 94 • C for 30 s, 56 • C at 1 min, 72 • C at 1 min, and a final extension at 72 • C for 5 min. To verify the pre-selective amplification success, 7 µL of the PCR products was electrophoresed on 2% agarose gel stained with GelRed at 100 V for 1 h. The pre-selective amplified PCR products were diluted 50X with ddH2O. Selective amplifications using four primers combination were performed using three labeled EcoRI primers (two primers labeled with FAM and one with PET (Table 3)) and four unlabeled Tru1I (Table 3). Selective primers were similar to the pre-selective ones with the addition of two nucleotides randomly chosen at the 3 end of the sequences ( Table 3). The selective PCR contained, for each reaction, 4 µL of the diluted pre-selective PCR product, 1X of Dream Taq buffer, 0.5 mM of MgCl 2 , 0.2 mM of each dNTPs, 0.2 µM of each selective primer, and 1 unit of DreamTaq polymerase (Thermo Fisher Scientific, Waltham, MA, USA) in a final volume of 25 µL. The PCR protocol was initial denaturation at 94 • C for 1 min, 10 cycles at 94 • C for 1 min, annealing at 65 • C to 56 • C (touchdown of 1 • C per cycle) for 1 min, 72 • C for 1 min 30 s, followed by 23 cycles at 94 • C for 30 s, 56 • C for 30 s, 72 • C for 1 min, and a final extension at 72 • C for 3 min. To check the selective PCR reactions success, 7 µL of the PCR products was electrophoresed on 2% agarose and gel stained with GelRed at 100 V for 1 h.

AFLP Scoring
AFLP products along with GS-500 LIZ size standard were run on a capillary sequencer (Applied Biosystems ® 3730XL, San Francisco, CA, USA) at GENTYANE Platform-INRA (Clermont-Ferrand, France). AFLP migration profiles were analyzed using GeneMapper v.5 (Thermo Fisher Scientific). The default AFLP settings of GeneMapper v.5 software for allele parameters, peak detection algorithm, peak quality, and quality flags were used. Fragment sizes between 50 bp and 500 bp were taken into consideration in the analysis. The allele's bins were scored as "0" if the peak height was ≤50, "1" if the peak height was ≥100, and by a "check" for the bins with the peak heights between 50 and 100. Manual analysis of the bin was performed based on the control individuals. Each peak that was present in the negative controls was discarded from the following analysis. Secondly, for each allele, genotyping error based on repetition type 1 (repetition from same DNA and independent restriction digests) and repetition type 2 (repetition from the same sample and different DNA extractions) was estimated using the Bonin error rate [28]. Loci that showed an error rate less than 10% were taken into consideration for further analyses. All loci that showed an error rate more than 10% were discarded. Then, peaks patterns that were scored as "check" were visually examined in GeneMapper and manually edited. All loci displaying singletons were also removed from the analyses.

AFLP Data Analysis
The proportion of polymorphic fragments (PLPs) for each primer combination was determined using AFLP-SURVEY 1.0 software [29]. Pairwise dissimilarity indices between individuals were calculated according to Nei and Li [30]. A neighbor joining (NJ) tree was created using PAUP v.4 based on the matrix of individuals dissimilarities. A heuristic search for the best NJ tree was conducted under the optimality criterion distance (minimum evolution (ME)) with 1000 bootstrap replicates. The starting seed was automatically generated. All characters had equal weight, and the branch-swapping algorithm followed tree bisection-reconnection (TBR). The builtNJ tree was unrooted and visualized using Interactive Tree Of Life (iTOL) v5 [31].
Because the occurrence of current or past events of interspecific hybridization has recently been documented between J. sabina and J. thurifera [23,24] or has involved an ancestral lineage of these species [22,32,33], we decided to also look at signals of genetic admixture between the species within the Sabina section represented in our sample. To achieve this goal, a Bayesian clustering approach was carried out on the basis of the AFLP fragments on the whole dataset using TESS software v. 2.3.1 [34]. This procedure allowed us to infer the number (K) of unobserved ancestral genetic clusters that best explained the sample genetic structure. When the model with admixture was chosen, the proportion of individual genomes originating from each potential ancestral genetic clusters (individual ancestry coefficients (q)) was estimated. In this case, the model took into account the possibility of historical events of genetic admixture of clusters. The statistical method implemented in TESS is spatially explicit in the sense that the admixture model assumes individual genetic relatedness are spatially auto correlated in a continuous manner (isolation by distance).
A first set of 140 runs (without admixture and correlated allele frequencies model) was carried out, varying K from 2 to 15 (10 runs per value of K) in order to estimate the best K values using the deviance information criterion (DIC) as a statistical variable for model choice ( Figure S1). MCMCs were carried out using a burn-in period of 70,000, followed by 120,000 iterations. The best values of K were chosen as those for which the plate with the smallest deviance information criterion (DIC) values was reached. Fifty new runs were performed for each of those K values using the admixture (CAR) model.
In order to identify potential different solutions for each K value (due to potential multimodality in the posterior distributions of the q values) and to manage the label switching among replicates in the TESS outputs, CLUMPP v1.1.2 was used with the Greedy option [35]. The probability of each solution was estimated based on their frequency among Life 2023, 13, 1479 7 of 17 50 runs. Plots representing the individual genomes ancestries to the inferred ancestral genetic clusters were created using the online STRUCTURE PLOT v.2.0 application. Because Bayesian clustering analyses on genotyping data always generate a background of spurious genome admixture in each individual, and because this background was mostly less than 10% in most individuals belonging to diploid species, individual ancestry coefficients less than 10% were not considered as an indication of interspecific genome admixture.

Polymorphism of AFLP Markers between Juniperus Species
In total, the 4 primer combinations generated 1022 loci used for the analyses. The proportion of polymorphic loci (PLPs) differed between species and primer combinations ( Table 4). The PLP ranged from 38.4% (for J. polycarpos) to 68.5% (for J. chinensis) in the primer combination EcoRI-1/Tru1I-1. In the primer combination EcoRI-1/Tru1I-3, the PLP ranged from 36.2% (for J. thurifera) to 58.5% (for J. chinensis). The EcoRI-2/Tru1I-4 primer combination showed the highest PLP among all the primer combinations, with 73.1% (for J. chinensis); the lower PLP in this primer combination was 43.8% (for J. thurifera). The primer combination EcoRI-4/Tru1I-2 showed a PLP ranging from 31.2% (for J. turcomanica) to 59.9% (in J. chinensis). Interestingly, J. chinensis showed the highest PLP for all primer combinations despite the low sample size analyzed (only nine individuals). Table 4. Proportion (PLP) of polymorphic loci (%) generated by AFLP from the four primer combinations of the studied Juniperus species. The red asterisk (*) represents the highest PLP in a primer combination, and the black asterisk (*) represents the lowest PLP.

Species
Primers EcoRI-1/Tru1I-1 This observation deserves to be explored on the basis of a more representative sampling scheme.

Phylogenetic Relationships among Juniperus Taxa
The NJ phylogenetic tree was supported with high bootstraps on all nodes and provided a similar pattern to that of Bayesian clustering. This showed that all individuals of the same species clustered together except in three cases ( Figure 1). Firstly, two individuals of J. excelsa from the Mrebbine population in north Lebanon belonged to the same evolutionary lineage of all J. polycarpos individuals. The second case concerned all J. turcomanica (=J. polycarpos var. turcomanica) individuals that were intermixed with J. polycarpos. The third case was observed for two individuals of J. polycarpos from the population Wadi El Njass in Lebanon collected in this study. These individuals clustered closely to the J. seravschanica lineage (Figure 1). polycarpos. The third case was observed for two individuals of J. polycarpos from the population Wadi El Njass in Lebanon collected in this study. These individuals clustered closely to the J. seravschanica lineage (Figure 1). By taking J. phoenicea as an outgroup, the NJ tree revealed that all studied taxa belong to two major lineages. The first one contains all individuals of J. chinensis, J. foetidissima, J. Sabina, and J. thurifera. The second lineage includes all individuals of J. excelsa, J. polycarpos, J. procera, J. seravschanica, and J. turcomanica. In both major lineages, each species is clearly separated from the others, except for J. polycarpos, J. turcomanica, and J. seravschanica. The latter species show unclear phylogenetic boundaries, despite the high number of AFLP markers studied. The current topology of the phylogenetic tree suggests that one or several ancestral lineages of J. thurifera, J. chinensis, and J. sabina could have contributed to the origin of the hexaploid J. foetidissima.

Admixture Patterns between Juniperus Taxa
The Bayesian clustering analysis provided K = 5, 6, and 7 ancestral genetic clusters as the best fits for the present data. However, when six and seven genetic clusters were considered, one and two "empty" clusters, respectively, were clearly identified (clusters with only very low or null coancestry values across the whole sample;). We therefore focused only on K = 5 as the best and more realistic number of ancestral genetic clusters. Three solutions for K = 5 (Figures 2 and S2) were identified on the basis of 50 runs. The most likely solution (Figure 2) held the highest probability (50%), followed by the less likely solutions A (36%) and B (14%) ( Figure S2). In addition, the most likely solution (Figure 2) was the only one that displayed five truly distinct genetic groups, because By taking J. phoenicea as an outgroup, the NJ tree revealed that all studied taxa belong to two major lineages. The first one contains all individuals of J. chinensis, J. foetidissima, J. Sabina, and J. thurifera. The second lineage includes all individuals of J. excelsa, J. polycarpos, J. procera, J. seravschanica, and J. turcomanica. In both major lineages, each species is clearly separated from the others, except for J. polycarpos, J. turcomanica, and J. seravschanica. The latter species show unclear phylogenetic boundaries, despite the high number of AFLP markers studied. The current topology of the phylogenetic tree suggests that one or several ancestral lineages of J. thurifera, J. chinensis, and J. sabina could have contributed to the origin of the hexaploid J. foetidissima.

Admixture Patterns between Juniperus Taxa
The Bayesian clustering analysis provided K = 5, 6, and 7 ancestral genetic clusters as the best fits for the present data. However, when six and seven genetic clusters were considered, one and two "empty" clusters, respectively, were clearly identified (clusters with only very low or null coancestry values across the whole sample;). We therefore focused only on K = 5 as the best and more realistic number of ancestral genetic clusters. Three solutions for K = 5 (Figure 2 and Figure S2) were identified on the basis of 50 runs. The most likely solution (Figure 2) held the highest probability (50%), followed by the less likely solutions A (36%) and B (14%) ( Figure S2). In addition, the most likely solution (Figure 2) was the only one that displayed five truly distinct genetic groups, because solutions A and B each displayed one "empty" genetic group, with the fifth one being ghosted as part of the genetic identity of J. chinensis (genetic cluster represented by the yellow in Figure S2). For the aforementioned reasons, we based our interpretation in this study mainly on the most likely solution presented in Figure 2.
Life 2023, 13, x FOR PEER REVIEW 9 of 17 solutions A and B each displayed one "empty" genetic group, with the fifth one being ghosted as part of the genetic identity of J. chinensis (genetic cluster represented by the yellow in Figure S2). For the aforementioned reasons, we based our interpretation in this study mainly on the most likely solution presented in Figure 2. The Bayesian clustering therefore revealed that J. phoenicea clearly belongs to a separate genetic cluster. Juniperus polycarpos, J. turcomanica, and J. seravschanica were overwhelmingly assigned to the same cluster. However, individuals from J. seravschanica showed an admixed genome composition between the J. polycarpos cluster and the same cluster as J. thurifera and J. foetidissima. Both varieties of J. sabina, i.e., J. sabina var. sabina (2n = 2x) and J. sabina var. balkanensis (2n = 4x), belong to the same genetic cluster with no admixture pattern, except for two individuals of J. sabina var. balkanensis, which presented a slight level of admixture with the genetic cluster of J. foetidissima and J. thurifera.
Juniperus excelsa and J. procera belong to the same genetic cluster. Remarkably, three individuals of J. excelsa were purely assigned to two different genetic clusters, two individuals to the J. polycarpos, cluster and the remaining one to the J. thurifera and J. foetidissima cluster. Additionally, a third J. excelsa individual displayed a strong signal of genome coancestry with the J. sabina cluster.
Juniperus chinensis showed a clear trend of ancestry with the J. sabina cluster, in accordance with their evolutionary proximity shown on Figure 1. However, it also displayed a complex mosaic pattern of coancestry with the genetic clusters of J. excelsa/J. procera and J. phoenicea.
The hexaploid J. foetidissima and the tetraploid J. thurifera were clearly assigned to the same genetic cluster. This was also very clear in solution B (only four genetic clusters, see above) and only partly confirmed in solution A ( Figure S2).
In the most likely solution (Figure 2), J. thurifera showed a slight signal of genetic admixture, with two individuals displaying genome coancestry (21% and 22%) with the same genetic cluster as J. sabina (a trend that was again more obvious in solution A but not in solution B). Similarly, these two genetic clusters were shown in the genetic ancestry of two J. foetidissima individuals from Greece (18% and 36%). In contrast, the population of J. foetidissima from Lebanon did not show any admixture signals in its ancestry.  [7] is/are represented by a circle colored in green, blue, and dark-red for diploid, tetraploid, and hexaploid levels, respectively.
The Bayesian clustering therefore revealed that J. phoenicea clearly belongs to a separate genetic cluster. Juniperus polycarpos, J. turcomanica, and J. seravschanica were overwhelmingly assigned to the same cluster. However, individuals from J. seravschanica showed an admixed genome composition between the J. polycarpos cluster and the same cluster as J. thurifera and J. foetidissima. Both varieties of J. sabina, i.e., J. sabina var. sabina (2n = 2x) and J. sabina var. balkanensis (2n = 4x), belong to the same genetic cluster with no admixture pattern, except for two individuals of J. sabina var. balkanensis, which presented a slight level of admixture with the genetic cluster of J. foetidissima and J. thurifera.
Juniperus excelsa and J. procera belong to the same genetic cluster. Remarkably, three individuals of J. excelsa were purely assigned to two different genetic clusters, two individuals to the J. polycarpos, cluster and the remaining one to the J. thurifera and J. foetidissima cluster. Additionally, a third J. excelsa individual displayed a strong signal of genome coancestry with the J. sabina cluster.
Juniperus chinensis showed a clear trend of ancestry with the J. sabina cluster, in accordance with their evolutionary proximity shown on Figure 1. However, it also displayed a complex mosaic pattern of coancestry with the genetic clusters of J. excelsa/J. procera and J. phoenicea.
The hexaploid J. foetidissima and the tetraploid J. thurifera were clearly assigned to the same genetic cluster. This was also very clear in solution B (only four genetic clusters, see above) and only partly confirmed in solution A ( Figure S2).
In the most likely solution (Figure 2), J. thurifera showed a slight signal of genetic admixture, with two individuals displaying genome coancestry (21% and 22%) with the same genetic cluster as J. sabina (a trend that was again more obvious in solution A but not in solution B). Similarly, these two genetic clusters were shown in the genetic ancestry of two J. foetidissima individuals from Greece (18% and 36%). In contrast, the population of J. foetidissima from Lebanon did not show any admixture signals in its ancestry.

Discussion
AFLP markers have been shown to display strong phylogenetic signals among closely related species [17][18][19][20][21]. Homoplasy or comigration of DNA fragments generated from different loci [29] is a potential challenge to the phylogenetic reconstruction based on AFLP markers. However, it was demonstrated to be a minor issue for studies dealing with closely related species [16], because homoplasy is expected to be less frequent. Additionally, the dominant nature of AFLP markers (the absence of a given amplified fragment being recessive) was also pointed out as a potential problem in assessing the genetic distance between a polyploid and each of its parents because the level of genotyping uncertainty increases with the ploidy level [36]. However, as noted by these authors, this should not inhibit the detection of hybridizations. Numerous studies have confirmed the power of AFLP markers in detecting interspecific hybridization (e.g., [23,24,[37][38][39]). Altogether, the 1022 loci produced in this study allowed a powerful discrimination between the studied species and the identification of several individuals that could be considered as progenies of interspecific hybridization.

Genetic Delimitation of Diploid Juniperus Taxa
Globally, the Bayesian clustering approach and phylogenetic analysis produce largely convergent results. Among the studied taxa, J. phoenicea, J. excelsa, J. procera, J. polycarpos, and J. turcomanica are exclusively diploid [7]. The results presented above showed that J. phoenicea is clearly an independent evolutionary lineage, consistent with available phylogenetic analysis inferring this species to the base of all ca. 60 species of the Sabina section [8].
Our data also showed that Juniperus excelsa and J. procera are genetically strongly related to each other, a result that is congruent with previous phylogenetic analyses based on ITS and four chloroplast genes [7,8]. Unexpectedly, two individuals of J. excelsa from the Mrebbine population in north Lebanon belong to the same genetic group as J. polycarpos in the Bayesian clustering and were located within the cluster of J. polycarpos in the NJ tree, indicating a possible misidentification of these individuals. Due to the morphological uniformity between the two species, J. polycarpos was discovered in Lebanon only few years ago based on genetic investigations [8,25,40,41]. Indeed, the separation between J. polycarpos and J. excelsa has been demonstrated solely based on microsatellite (SSR), nrDNA (ITS), and chloroplast sequences [40,41] and never on the basis of morphological traits. Therefore, along with the mentioned markers, the AFLP markers provided a very clear and strong separation between these two species and represent a future promising cost-effective method for a heuristic geographic screening of these taxa. On the other hand, the AFLP markers failed to differentiate between J. polycarpos and J. turcomanica (=J. polycarpos var. turcomanica). In the present study, both species belonged to the same cluster without a clear delimitation between them (Figures 1 and 2). Congruently, previous studies using a few nuclear (low copy genes and ITS) and four chloroplast sequences [41,42] showed no clear separation between these two taxa. Therefore, based on the current and aforementioned studies, we propose considering J. polycarpos and J. turcomanica as a species complex. Deep phylogenomic analyses must be conducted to obtained a clear understanding of this complex species, notably using NGS techniques such as GBS [43] and Hyb-Seq [44].

Potential Origin of the Tetraploid Junipers
Juniperus sabina, J. chinensis, and J. seravschanica have two cytotypes (2n = 2x and 4x) [7], whereas J. thurifera was the only studied species that is exclusively tetraploid. In this study, we analyzed two varieties of J. Sabina: the diploid J. sabina var. sabina (population: Austria, Alps) and the tetraploid J. sabina var. balkanensis (populations: Bosnia-Herzegovina and Croatia). Both varieties clustered in the same branch in the NJ tree and belonged to the same genetic group in the Bayesian clustering analysis. Only very few events of admixture were found for the tetraploid variety, and this result could be interpreted as a signal of shared ancestral polymorphism. Recently, J. sabina var. balkanensis has been described as a hybrid between the diploid variety of J. sabina and the tetraploid J. thurifera, since it holds the nuclear nrDNA (ITS) sequence of J. sabina var. sabina and the chloroplast haplotype (based on four genes) of J. thurifera [22,32,33]. These studies failed to uncover any evidence of a contribution of J. thurifera to the nuclear marker (ITS) of J. sabina var. balkanensis, in agreement with the AFLP in the present study despite the large number of loci. Therefore, it is likely that, after the ancient hybridization between a female J. sabina and a male J. thurifera, at least most of parental genetic nuclear material from J. thurifera was removed in the lineage of J. sabina var. balkanensis. This could have been achieved by several generations of backcross to the female J. Sabina, as described in one of the parsimonious hypothetical polyploidization pathways presented by Farhat et al. [12]. However, meiotic disorders may have helped in accelerating this process of elimination of one parental genetic material from the genome of the interspecific hybrid. Recently, interspecific hybridization was detected between the J. sabina (2n = 2x) and J. thurifera (2n = 4x) present in sympatric populations in France (Alps) and Spain (in two populations), giving rise to triploid and possibly tetraploid individuals [23,24]. Interestingly, male triploid hybrids examined in the French Alps population showed well-conformed and potentially fertile pollen [23], which might reflect the fertility of the triploid hybrids. Consequently, the presence of fertile triploid hybrids might engender a better establishment of the polyploid toward a more stable ploidy level. These latter findings provide more support for the hybrid origin of the tetraploid J. sabina var. balkanensis involving the parental species J. thurifera and J. sabina.
Juniperus chinensis showed a complex genetic pattern. However, interpretation of this pattern should be very cautious. It could be partly due to the fingerprint of shared ancestral polymorphism and/or to uncertainties in ancestry inferences because of the constrained and limited number of ancestral genetic cluster of the solutions retained after Bayesian clustering (five or four genetic clusters in this analysis). Nevertheless, our results strongly underline the strong coancestry of J. chinensis individuals with J. sabina individuals. In addition, J. chinensis was located in the NJ tree on the sister branch of J. sabina. This is compatible with a previous analysis showing that J. sabina and J. chinensis are phylogenetically closer to each other than the remaining taxa examined in our study [8]. Additionally, the complex admixed genetic pattern of J. chinensis found in this study might be the consequence of plausible introgression events from neighboring taxa. Juniperus chinensis is currently distributed in China and Japan, far from all the species analyzed in this study, except for J. sabina var. sabina, which has a wide distribution reaching China ( Figure S3) [8]. Therefore, if the observed complex genetic pattern in J. chinensis was produced by introgression events, it might implicate historically dated events during which the geographical distribution of the studied species was possibly closer than that observed currently.
The revealed pattern of genetic admixture of the genome of most J. seravschanica individuals between the genetic groups of J. polycarpos/J. turcomanica and J. foetidissima/J. thurifera suggests that this species is a hybrid. This result is in accordance with the phylogenetic closeness previously revealed between J. seravschanica and J. foetidissima [8]. Moreover, it has lately been suggested based on nrDNA (ITS) and four chloroplast regions that J. seravschanica obtained the chloroplasts from an ancestor of J. foetidissima/J. thurifera through an ancient chloroplast capture [45], meaning that hybridization has occurred in the past between these two lineages (J. seravschanica and J. foetidissima/J. thurifera). Interestingly, the studied populations of J. seravschanica were previously identified as tetraploid [7]. Therefore, the AFLP data suggest that tetraploid individuals of J. seravschanica issued from an interspecific hybridization event, and the inferred admixed status of the genome of some of J. seravschanica individuals may testify to their allopolyploid nature. Based on the admixed pattern revealed in this study and the chloroplast capture analysis of Adams [45], we propose that the J. polycarpos/J. turcomanica lineage has contributed maternal parent and the J. foetidissima/J. thurifera lineage has contributed the paternal parent of tetraploid J. seravschanica. In addition, the strong admixture level shown by a few J. seravschanica individuals (ca. 50% of the admixture) ( Figure 2) may also be the consequence of still ongoing genetic introgression. Indeed, the current geographical distributions of J. foetidissima, J. polycarpos, J. turcomanica, and J. seravschanica are currently overlapping. More in-depth analysis using large and genome-wide DNA sequence data should be conducted with a more extensive and more geographically representative sample to look further into these hypotheses. Especially, the sympatry zones of J. seravschanica with parental species candi-somal segregation may produce diverse genetic combination in the progenies of newly formed polyploid individuals with variable contribution of genomes from both parental species (e.g., [46]). Additionally, selection can then favor or eliminate certain combinations, or drift processes and the founder effect along pathways of dispersal can lead to the random fixation of certain genetic combinations. All those processes could have happened in the case of the evolution of J. foetidissima in the eastern part of the Mediterranean region.
( Figure 3, pathway I). A whole-genome duplication (WGD) in triploid hybrids allows us to restore cytogenetic stability and full fertility and is therefore selectively advantageous (Figure 3, pathway I). According to this pathway, the expected genomic contribution of J. sabina to the hexaploid genome of J. foetidissima would be one-third at the beginning of the process, which is in accordance with the admixed genome pattern of only a few individuals from Greece in our study. However, different mechanisms such as high rates of nonhomologous association, homologous recombination, homologous chromosome replacement, and bias in chromosomal segregation may produce diverse genetic combination in the progenies of newly formed polyploid individuals with variable contribution of genomes from both parental species (e.g., [46]). Additionally, selection can then favor or eliminate certain combinations, or drift processes and the founder effect along pathways of dispersal can lead to the random fixation of certain genetic combinations. All those processes could have happened in the case of the evolution of J. foetidissima in the eastern part of the Mediterranean region. Alternative pathways are illustrated in Figure 3 (pathways II and III), following possibly similar scenarios as discovered in the French Alps and Spain [23,24]. The triploid generated after the cross between ancestors of J. thurifera (n = 2x) and J. sabina (n = x) would, in turn, produce gametes with three ploidy levels: n = x, 2x, and 3x. It is very unlikely that the triploid gamete would interfere in this scenario with the ancestor of J. thurifera (2n = 4x) because it would have resulted in pentaploid progenies. Therefore, we discard this possibility. In the first case, the cross of the (n = x) gamete of the neo-triploid hybrid with (n = 2x) of J. thurifera ancestor would produce another triploid hybrid that would duplicate its genome through a WGD reaching the hexaploid level ( Figure 3, pathway II). In the second case, the diploid gamete (n = 2x) of the neo-triploid hybrid might unify with the J. thurifera ancestor (n = 2x), giving rise to tetraploid level progenies, as found in the French Alps ( Figure 3, pathway III). A further hybridization between the triploid and tetraploid progenies (2n = 3x and 2n = 4x) may give rise to triploid hybrids (among other ploidy Alternative pathways are illustrated in Figure 3 (pathways II and III), following possibly similar scenarios as discovered in the French Alps and Spain [23,24]. The triploid generated after the cross between ancestors of J. thurifera (n = 2x) and J. sabina (n = x) would, in turn, produce gametes with three ploidy levels: n = x, 2x, and 3x. It is very unlikely that the triploid gamete would interfere in this scenario with the ancestor of J. thurifera (2n = 4x) because it would have resulted in pentaploid progenies. Therefore, we discard this possibility. In the first case, the cross of the (n = x) gamete of the neo-triploid hybrid with (n = 2x) of J. thurifera ancestor would produce another triploid hybrid that would duplicate its genome through a WGD reaching the hexaploid level ( Figure 3, pathway II). In the second case, the diploid gamete (n = 2x) of the neo-triploid hybrid might unify with the J. thurifera ancestor (n = 2x), giving rise to tetraploid level progenies, as found in the French Alps (Figure 3, pathway III). A further hybridization between the triploid and tetraploid progenies (2n = 3x and 2n = 4x) may give rise to triploid hybrids (among other ploidy levels), followed by a WGD reaching the hexaploid level. In both pathways II and III, the original contribution of the genome of J. sabina to the genome of J. foetidissima may widely vary from zero to one-third.
On the other hand, if the hybridization between the ancestors of J. thurifera and J. sabina happened before the WGD in the J. thurifera lineage, the only scenario to reach a hexaploid level after hybridization would be through the contribution of "unreduced gametes" in one of the parents. This scenario would generate triploid hybrids followed by WGD. However, this hypothesis is not very plausible because the polyploidy event in J. thurifera was already estimated to be very old (paleo-tetraploid) based on cytological evidence [11].
All of the hypotheses discussed above suggest the diploid J. Sabina ancestor as one of the contributors to the J. foetidissima genome. The alternative possible hypothesis would