Asymmetric Reproductive Barriers and Gene Flow Promote the Rise of a Stable Hybrid Zone in the Mediterranean High Mountain

Hybrid zones have the potential to shed light on evolutionary processes driving adaptation and speciation. Secondary contact hybrid zones are particularly powerful natural systems for studying the interaction between divergent genomes to understand the mode and rate at which reproductive isolation accumulates during speciation. We have studied a total of 720 plants belonging to five populations from two Erysimum (Brassicaceae) species presenting a contact zone in the Sierra Nevada mountains (SE Spain). The plants were phenotyped in 2007 and 2017, and most of them were genotyped the first year using 10 microsatellite markers. Plants coming from natural populations were grown in a common garden to evaluate the reproductive barriers between both species by means of controlled crosses. All the plants used for the field and greenhouse study were characterized by measuring traits related to plant size and flower size. We estimated the genetic molecular variances, the genetic differentiation, and the genetic structure by means of the F-statistic and Bayesian inference. We also estimated the amount of recent gene flow between populations. We found a narrow unimodal hybrid zone where the hybrid genotypes appear to have been maintained by significant levels of a unidirectional gene flow coming from parental populations and from weak reproductive isolation between them. Hybrid plants exhibited intermediate or vigorous phenotypes depending on the analyzed trait. The phenotypic differences between the hybrid and the parental plants were highly coherent between the field and controlled cross experiments and through time. The highly coherent results obtained by combining field, experimental, and genetic data demonstrate the existence of a stable and narrow unimodal hybrid zone between Erysimum mediohispanicum and Erysimum nevadense at the high elevation of the Sierra Nevada mountains.

Introduction

Evaluating the mechanisms promoting the rise of hybrid zones helps us to understand the nature and dynamics of these interesting evolutionary scenarios. Hybrid zones have been pointed to as one of the most insightful places to study the evolutionary interactions between divergent but related taxa (Barton and Hewitt, 1985; Harrison, 1993). Important effects of hybrid zones were well documented on, for example, morphological traits (Tastard et al., 2012, Keller et al., 2021), behavior (Good et al., 2000), selective regimes (Cruz et al., 2001; Knief et al., 2019), and biodiversity (Whitham et al., 1994, 1999). The hybrids can show a modified phenotype promoting their evolvability and, by altering the reproductive barriers between the parental species, modify the evolvability of these lasts (Parsons et al., 2011). However, an integrative study of the mechanisms, patterns, and consequences of hybridization is not straightforward.

Even before the first conceptual models to explain hybrid zones were proposed, gene transfer between the neighboring species had been recognized as an important evolutionary process (Dobzhansky, 1940; Mayr, 1942; Anderson, 1949). The main classifications of hybrid zones proposed relate to the distribution of the hybrid and the parental phenotypes/genotypes (Harrison and Bogdanowicz, 1997; Jiggins and Mallet, 2000). Therefore, hybrid zones can be classified as unimodal, bimodal, trimodal, and asymmetric (Pickup et al., 2019). Unimodal hybrid zones are characterized by a higher frequency of individuals with intermediate genotypes between both parental species (Jiggins and Mallet, 2000). Hybridization and admixture are predominant in this kind of hybrid zones, so weak reproductive barriers are expected in them (Figure 1). A special case of the unimodal hybrid zones appears when the hybrid zone forms a single panmictic population exhibiting individuals with variable degrees of genetic similarity to the parental forms, termed as “hybrid swarms” (Harrison and Bogdanowicz, 1997; Jiggins and Mallet, 2000). Hybrid swarms could present individuals coming from first or second hybrid generations or from different levels of back-crossing. In a bimodal hybrid zone, individuals belonging to the parental species co-occur with a very low frequency of their hybrids (Cruzan and Arnold, 1993; McMillan et al., 1997; Jiggins and Mallet, 2000; Vedenina and Helversen, 2003). This kind of hybrid zone appears when strong reproductive isolation barriers exist between the species previously (Harrison and Bogdanowicz, 1997; Jiggins and Mallet, 2000; Coyne and Orr, 2004; Figure 1). However, when reproductive isolation between taxa is only partial, a stable hybrid zone would arise, generating a trimodal hybrid zone that depends, on the one hand, on the differences in fitness between the parentals and the hybrids and, on the other hand, on the migration rates from the parental zones (Key, 1968; Barton and Hewitt, 1985; Jiggins and Mallet, 2000; Figure 1). Finally, asymmetric hybrid zones occur when strong asymmetries in gene flow and barrier strength promote hybridization and introgression in one direction, such that gene flow occurs only from one species to another (Pickup et al., 2019; Figure 1).

FIGURE 1

Figure 1. (A) Different models of hybrid zones based on the different patterns of the genetic mixture between divergent taxa. The abscissa axis shows individuals in the hybrid zone and the ordinate axis shows the ancestry proportion exhibited by individuals. (B) Frequency of each range of ancestry proportion in the population and the gene flow expected in each situation. Blue, green, and black arrows represent the gene flow promoted by interspecific crosses, backcrosses, and crosses between hybrids, respectively. The persistence of gene flow is indicated by solid lines; the absence of gene flow is indicated by broken lines. (C) Patterns of the phenotypic differentiation expected when differences between hybrids and parental plants are found. H1 and H2 are the hybrids obtained by adding pollen on parental plants 1 and 2, respectively. Yellow and red dots and bars correspond to the parental species, while orange represents hybrid values.

Hybridization is much more frequent and is evolutionarily relevant in plants than in animals, but surprisingly, in the past, it has been more explored in the latter (Ellstrand et al., 1996; Dowling and Secor, 1997). The number of plant hybrid zones studied using an integrative approach, including genetic analyses and the nature of reproductive isolation, is limited (Abbott, 2017). However, interesting examples of contact zones and hybrid populations have been described in plants, for example, those described in Phlox (Levin, 1967), Iris (Cruzan and Arnold, 1993; Young, 1996; Emms and Arnold, 1997), Helianthus (Rieseberg et al., 1998), or Armeria (Fuertes Aguilar et al., 1999) and, more recently, those described in Ipomopsis (Campbell and Aldridge, 2006), Silene (Minder et al., 2007), Narcissus (Marques et al., 2010), or Primula (Keller et al., 2021), among others. These studies have demonstrated how the rise of hybrids can not only generate new phenotypic values in the contact zone but also modify the parental phenotype through introgression (Owens et al., 2016; Nieto Feliner et al., 2019) or by affecting the selective trajectories in the contact zone and the surrounding areas (Wielstra et al., 2017; Wielstra, 2019). Therefore, it is interesting to explore the phenotypic differences between the parental plants and hybrids not only in natural populations but also in controlled conditions where we could establish the origin and the level of admixture of the hybrid plants. In this sense, the hybrids can show different patterns of phenotypic differentiation. The hybrid phenotype could exhibit higher, lower, or intermediate phenotypic values compared with the parental plants (Figure 1C). The hybrid can also show an asymmetric pattern when the hybrids on one parental direction (that is, the parental direction acting as the mother) exhibit significant differences with the hybrid in the other sense (Figure 1C). Exploring this question through time in natural populations would allow us to know more about the stability of hybrids and the hybrid zones. However, the longer the life cycle of the organism being experimented on, the more difficult exploring becomes.

Incomplete reproductive barriers (or absence of them) are necessary to the occurrence of hybrids in a contact zone between different species. So, the evaluation of the strength of reproductive barrier components between the hybridizing species sheds light on the mechanisms underlying the rise of hybrid zones. Again, this objective is not always easily reachable, especially when the study systems demand several to many years for blooming. Reproductive barriers are usually estimated by the use of one or a few components of fitness. However, quantifying the different fitness components through the entire life cycle of plants will allow for a more realistic estimation of the intensity of these barriers and their cumulative effects on gene flow between the parental species and hybrids (Baack et al., 2015).

In this study, we examined two plant species belonging to the genus Erysimum L., Erysimum mediohispanicum, and Erysimum nevadense, and the hybrids produced between them in the Sierra Nevada (Southeast Spain). Genus Erysimum L. (Brassicaceae) is mainly distributed in the northern hemisphere (Polatschek, 1986), presenting an important diversification center in the western Mediterranean region (Greuter et al., 1986). The genus presents a complex evolutionary history due to events of interspecific hybridization and polyploidization, producing a highly diversified genus (more than 200 species) enriched in species complexes and cryptic species (Clot, 1992; Ancev, 2006; Marhold and Lihová, 2006; Turner, 2006; Couvreur et al., 2010, Abdelaziz et al., 2011; Al-Shebaz, 2012). In the Sierra Nevada, E. mediohispanicum and E. nevadense show contrasting distributions, inhabiting the lowland and the top of the mountains, respectively (Blanca et al., 2009). However, these species make what could be a secondary contact (Abdelaziz et al., 2014) at mid-elevation.

In the present study, we used an integrative approach by combining fieldwork, lab experiments, and molecular data to evaluate the nature and the strength of barriers to hybridization, to characterize the hybrid zone between E. mediohispanicum and E. nevadense in the Sierra Nevada National Park, and to explore the consequences of this geographic contact on plant phenotypes. The main goals of this study are (1) to explore the genetic and phenotypic differences between both species and their hybrids in the contact zone; (2) to explore the gene flow patterns promoting hybridization in this area; (3) to analyze the stability of the contact zone through time; and (4) to analyze the reproductive barriers promoting or constraining the rise of hybrids between E. mediohispanicum and E. nevadense.

Materials and Methods

Study System

Erysimum mediohispanicum Polatschek is an endemic species of the Iberian Peninsula, where it is distributed in two widespread and almost disconnected regions, located in the North and Southeast of the peninsula. Its life cycle varies among individuals and populations, usually being monocarpic. In the Sierra Nevada (Southeast of Spain), E. mediohispanicum is composed of diploid populations (Muñoz-Pajares et al., 2018) of facultative biennial, spending 2–3 years growing at an elevation of up to 2,200 m like a vegetative rosette on calcareous soils (Muñoz-Pajares et al., 2020). After that period, plants display flowers on one to three long stalks, ranging in number from only a few to several hundred (Gómez, 2003). These flowers are visited by a highly diverse assemblage of insects (Gómez et al., 2007; Gomez et al., 2014).

Erysimum nevadense Reut. is also a diploid species (author's unpublished data) and mostly polycarpic herb endemic to the peaks of the Sierra Nevada mountains, where it grows on siliceous soils at an elevation from 2,300 to 2,700-m (Blanca et al., 2009). E. nevadense spends 2–3 years like a rosette before displaying anything from a few to several hundred flowers on various short stalks. It is also a pollination-generalist plant, but it does not present a pollinator assemblage as diverse as that of E. mediohispanicum, probably due to the harsh conditions of its habitat (Gómez et al., 2007; Mohamed, 2013).

Both species come into contact along a narrow area at an elevation of approximately 2,200 m. In this contact zone, populations from each species are located at a mere hundred meters apart, generating an area where the phenotypic traits defining each species are difficult to identify. In 2007, we established a line transect of 3.1 km between these species in the north face of the Sierra Nevada encompassing five populations (Table 1), two E. mediohispanicum populations (Em25 and Em17), two E. nevadense populations (En11 and En10), and a putative hybrid population (H01), where both species contact. After 10 years, in 2017, we phenotypically resampled the hybrid population (H01) and the two parental populations more genetically related to the hybrid one (Em25 and En10).

TABLE 1

Table 1. Geographic information and sampling effort in the five populations included in the study (ordered by elevation).

Plant Phenotyping

For about a month from the end of May in 2007, at each of these populations, we phenotyped 90 plants by measuring the following phenotypic traits: (1) number of flowers and (2) number of stalks: counting the total number of flowers and stalks produced by each plant, respectively; (3) stalk height: the height of the tallest stalk of the plant from the ground to the top of the stalk at the end of the flowering period, when no more flowers are expected on the stalk, using a measuring tape (±0.5 mm error); (4) stalk diameter: the basal diameter of the tallest stalk; (5) corolla diameter: the distance between the edge of the two opposite petals; (6) corolla tube width: the diameter of the corolla tube aperture as the distance between the bases of two opposite petals; and (7) corolla tube length: the distance between the corolla tube aperture and the base of the sepals. The stalk diameter and the traits related to the flower size (4–7, above) were quantified in millimeters using a digital caliper (±0.01 mm. error). On approximately the same dates in 2017, all these phenotypic traits were remeasured in the same number of plants in the populations Em25, H01, and En10 (Table 1).

For each population and sampled year, we calculated the product–moment correlations, the covariances between each couple of the phenotyped traits, and the variance per trait (Supplementary Tables 1, 2). Kruskall–Wallis tests and one-way ANOVAs were used to compare the trait distributions and means of the populations for each phenotypic character, respectively. These analyses were performed using the package stats in R (R Development Core Team, 2014).

DNA Isolation and Genotyping

In 2007, we collected fresh tissues from all the phenotyped plants at the putative hybrid area and also from randomly selected 30 plants out of the 90 phenotyped plants from each of the four parental populations in the transect (in total, 210 plants). The tissue was stored in silica gel for the subsequent DNA isolation, using the GenElute Plant Genomic DNA Miniprep kit (Sigma-Aldrich, St. Louis, MO, USA). The isolated DNA was used for the individual genotyping, carried out with 10 microsatellite markers previously described by Muñoz-Pajares et al. (2011). We performed PCR in a 15-μL reaction mixture containing 0.17 ng/μL of the template genomic DNA, 1X buffer (ref. M0273S, New England BioLabs), 0.16 mM each of dNTP (Sigma-Aldrich), 0.33 μM each of forward (fluorescently tagged) and reverse primers, and 0.02 U/μL Taq polymerase (ref. M0273S, New England Biolabs). PCR was conducted in a Gradient Master Cycler Pro S (Eppendorf, Hamburg, Germany) with an initial step of 30 s of denaturation at 94°C followed by 35 cycles at 94°C for 15 s, annealing temperatures for each microsatellite marker described in Muñoz-Pajares et al. (2011) for 30 s, extension at 72°C for 30 s, and a final extension step at 72°C for 3 min. PCR products were diluted at a ratio of 1:15 and sent to MACROGEN (Geumchun-gu, Seoul, South Korea; http://www.macrogen.com) for capillary electrophoresis using 400HD ROX as standard. Alleles were called using the Peak Scanner Software version 1.0 (Applied Biosystems).

Genetic Data Analyses

Once we established the multilocus genotypes, we excluded the individuals presenting levels of missing data higher than 30% from the genetic analyses. This included 3 individuals from the En11 population and 10 individuals from H01 (Table 2). To genetically characterize the populations, we estimated the following parameters: (a) number of non-redundant multilocus genotypes (N_G), as the number of genotypes showing at least one different allele, excluding missing data; (b) mean number of alleles per locus (n_a); (c) observed heterozygosity (Ho), as the actual frequency of heterozygous individuals in the sample. We estimated the mean of the individual heterozygosities per population using the ratio between the number of heterozygote loci and the number of successfully genotyped loci; (d) gene diversity (Hs), as the expected proportion of heterozygous individuals assuming the Hardy–Weinberg equilibrium. Gene diversity was calculated using the Nei (1987) estimator; (e) mean allelic richness per locus (Rs), estimated as the probability of sampling the allele i at least once among the 2n genes of a sample, being independent of sample size; (f) the mean private allelic richness (R_P), estimated as the mean number of singular alleles per locus presented at each population. All the previous parameters were calculated using the package hierfstat v. 0.04–6 (Goudet, 2005) or using scripts developed by ourselves, both in R (R Development Core Team, 2014). (g) Inbreeding coefficient (F_IS), which provides information about the Hardy–Weinberg equilibrium departures due to either excess or defect of heterozygotes. We estimated F_IS by Bayesian inference using BayesAss v3.0 (Wilson and Rannala, 2003) for each population and overall for the two studied species. Analysis lasted for 10 million MCMC iterations, sampling every 1,000 generations and optimizing the mixing parameter for allele frequencies and for inbreeding coefficients. After that, we removed the first 10% of total iterations and checked the trace files with the program Tracer v1.4 (Rambaut and Drummond, 2007) to determine the convergence of the independent Bayesian MCMC runs.

TABLE 2

Table 2. Genetic diversity parameters and effective population sizes for Erysimum mediohispanicum (Em), Erysimum nevadense (En), and the hybrid populations (total N = 197).

Microsatellite-based genetic differentiation among groups of populations belonging to the same species, among populations within groups, and among individuals within populations was estimated using a hierarchical analysis of molecular variance (AMOVA), as implemented in Arlequin (Excoffier and Lischer, 2010), using 1,000 permutations to test the significance. We performed this analysis two times: the first one excluding the hybrid population and the second one including the hybrid population in its own group. In addition, pairwise comparisons for genetic differentiation between populations were computed using the package hierfstat in R (R Development Core Team, 2014) by the computation of FST statistics (Weir and Cockerham, 1984) and D_ST (Jost, 2008). F_ST significance was calculated from 1,000 permutations, while the DST was estimated as the harmonic mean across loci.

The genetic relationship among genotypes was inferred by Bayesian means using the model-based clustering algorithm, as implemented in Structure v.2.2 (Pritchard et al., 2000; Falush et al., 2003). The number of multilocus genotype clusters (K) was analyzed using diploid setting and using admixture and prior information as ancestry models and correlations as the allele frequency model using prior information. We performed simulations with 10 replicates for each K value, ranging from K = 1 to K = 6. Each run consisted of 50,000 Markov Chain Monte Carlo (MCMC) steps after 20,000 burn-in steps. To detect the optimum value of K, we applied the Evanno method (Evanno et al., 2005), as implemented in the Structure Harvester website (Earl and von Holdt, 2012).

Finally, we estimated the gene flow rates among the studied populations and their significance by means of Bayesian inference using BayesAss v3.0 (Wilson and Rannala, 2003). Analysis lasted for 10 million MCMC iterations, with a sampling frequency of every 1,000 generations, optimizing the mixing parameter for allele frequencies and for inbreeding coefficients. After that, we removed the first 10% of the total iterations and checked the trace files with the program Tracer v1.4 (Rambaut and Drummond, 2007) to determine the convergence of the independent Bayesian MCMC runs.

Reproductive Barrier Evaluation

In September 2009, we collected 120 individual plants from the four parental populations (30 plants per population of E. mediohispanicum and E. nevadense; Table 1). These individuals were transplanted to individual pots (11 cm x 11 cm x 15 cm) using the same soil in which they were growing and moved to a common garden in the University of Granada (an elevation of ~700-m). By May 2010, 44 individuals had survived (30 E. mediohispanicum plants and 14 E. nevadense plants; Table 1). Before the plants started blooming, they were moved to a greenhouse to exclude them from pollinators. During flowering, they were phenotyped following the abovementioned methodology used for the field plants, and several additional flowers per individual were subjected to two different treatments: (a) outcrossing (OC), in which some flowers were emasculated before opening and were pollinated with pollens from a different conspecific individual; and (b) interspecific hybridization (HC), in which some flowers were emasculated before opening and were pollinated with pollens from an individual from the other species. In total, 232 flowers were used in the experiment, with a mean of 5.5 ± 3 flowers per plant (3.2 ± 1.9 for OC and 2.8 ± 1.7 for HC per plant).

Once the blooming period was over, we recorded the number of flowers per plant and the treatment setting as ripe fruits or aborting without producing any fruits. The total number of ovules, unfertilized ovules, aborted seeds, and ripe seeds produced per ripe fruit were recorded in the lab using magnifying glasses. A total of 3,139 seeds were harvested at the end of the experiment. Subsequently, when possible, 15 seeds per plant and treatment were taken at random and sown randomly in a greenhouse. Their germination was recorded two times a week for the first month, and their survival was recorded every month during the next 10 months.

The following predispersal and postdispersal components of the plant reproductive output were quantified for each treatment and plant: (a) fruit set, the proportion of labeled flowers setting fruit; (b) seed production, the number of seeds produced per ovule in a given fruit; (c) seedling emergence, the proportion of sown seeds germinating and emerging as seedlings; and (d) seedling survival, the proportion of seedlings surviving until the end of the experiment. Afterward, we calculated the cumulative pre-dispersal fitness (W_PRE), as fruit set x seed production, and the cumulative total fitness (W_TOT), as fruit set x seed production x seedling emergence x seedling survival. The F1 generation resulting from OC and HC was grown until blooming (in the spring of 2011) when it was also phenotyped as we did with the natural populations.

Hybrid inviability (HI) was calculated per population by comparing the fitness between the intraspecific and interspecific crosses using the Ågren and Schemske (1993) approach as

$\begin{array}{ll}HI=1-w_{\text{h}}/w_{\text{o}};w_{\text{h}}<w_{\text{o}}& (1)\end{array}$

$\begin{array}{ll}HI=w_{\text{h}}/w_{\text{o}}-1;w_{\text{h}}>w_{\text{o}}& (2)\end{array}$

where w_o is the fitness of the intraspecific outcrossing treatments and w_h is the fitness of the hybrid crosses. In all cases, we used both predispersal and total fitness. Using this approach, the values range between −1 and +1. Positive values indicate that the hybrid crosses have lower fitness than intraspecific crosses (occurrence of hybrid inviability and hence the rise of reproductive barriers). The significant values of HI were calculated by computing 95% CIs by means of bootstrapping with 1,000 permutations, using package boot in R (Canty and Ripley, 2017).

Results

Phenotypic Differentiation

Using one-way ANOVAs, we found significant differences at the among-population level in the 2007 sample for all the measured phenotypic traits, except for the corolla tube width (Figure 2A). The hybrid population showed the highest values for the number of flowering stalks and the number of flowers, both traits related to plant size. However, it showed intermediate values for the rest of the traits, including characters related to plant size (stalk diameter and plant height) and to flower size (Figure 2A). The same pattern was found again when we resampled three out of the original populations in 2017 (Em25, H01, and En10; Table 1 and Figure 2C).

FIGURE 2

Figure 2. Values of phenotypic traits at each population in 2007 (A), 2017 (C), and the hybrid and parental plants obtained by controlled cross experiments (B). Yellow dots denote plants or populations of Erysimum mediohispanicum, orange dots denote the hybrid plants or population, and red dots denote plants or populations of Erysimum nevadense. F-ratios refer to one-way ANOVAs. Letters indicate the groups of plants where the differences are significant, according to Tukey's HSD comparison.

The analysis of phenotypic traits performed on the offspring (F1 generation) of the controlled crosses confirmed these trends but also showed different patterns depending on the direction of gene flow producing the hybrid plant. The plants resulting from the hybrid crosses performed on E. mediohispanicum, i.e., E. mediohispanicum acting as the female, showed significantly higher values than the parental species for the number of flowering stalks and the number of flowers, while the hybrid plants resulting from crosses on E. nevadense showed lower values for these same traits (Figure 2B). For the rest of the phenotypic traits, the hybrid plants presented intermediate values, not significantly different from the values shown by the parentals for stalk diameter, flower diameter, and flower tube width (Figure 2B).

We found significant correlations in every population between the phenotypic traits defining the plant size (number of stalks, stalk diameter, plant height, and the number of flowers) and between the phenotypic traits defining the flower size (corolla diameter, corolla tube length, and corolla tube width) (Supplementary Table 1). In some cases, we also found significant positive correlations between the plant-size traits and the flower-size traits. This means that larger plants can also produce larger flowers, but larger plants do not present smaller flowers in any case. This pattern was found for the plants measured in the natural populations in 2007 (Supplementary Table 1) and 2017 (Supplementary Table 2), as well as in the plants obtained by the controlled crosses. The plants obtained by the intraspecific crosses presented similar number and levels of correlations to the plants phenotyped in natural populations. The experimental within-species cross plants showed similar patterns to those found in plants produced from hybrid crosses, but E. nevadense hybrids, where the number of traits significantly correlated, decreased (Supplementary Table 3).

Genetic Diversity and Differentiation

A total of 197 plants were successfully genotyped from the five sampled populations in 2007 (Table 2). A total of 196 multilocus genotypes were identified from them, finding no identical multilocus genotypes among populations, but only one repeated multilocus genotype at population En11 (Table 2).

The mean number of alleles per locus (n_a) ranged between 6.0 and 9.9, corresponding to En11 and H01, respectively, and E. mediohispanicum presented slightly higher values than E. nevadense (Table 2). The observed H_O presented values from 0.50 to 0.63, corresponding to these extremes to En11 and Em25, respectively, and again E. mediohispanicum values being higher. Moreover, the mean H_S presented the same pattern shown by H_O but the values of H_S were higher than Ho in all the considered populations and species (Table 2). The values of the mean allelic richness (R_S) ranged from the minimum value (5.87) at En11 to the maximum (7.54) at H01, presenting higher values for E. mediohispanicum populations (Table 2). The pattern exhibited by the mean private allelic richness (R_P) was different, being higher for E. nevadense and presenting their minimum value at Em17 (0.10) and the maximum value at H01 (1.20). Finally, the inbreeding coefficient (F_IS) exhibited its minimum value at the hybrid population (0.139) and the maximum values at both extremes of the transect (Em25 = 0.363 and En10 = 0.360), although the SEs associated with these maximum values were also higher (Table 2).

Hierarchical analysis of molecular variance (AMOVA) indicated that population groups (species) were not genetically differentiated for these markers (Table 3). Moreover, this also occurred when H01 was included as a separate group. In contrast, the among-population-within–groups level (both including and excluding H01 in the analysis) was significantly differentiated, although it accounted for only low amounts of genetic variance (Table 3). The genetic variance at this level decreased when the hybrid population was included (Table 3). The within-population level accounted for almost all genetic variance in both cases (Table 3).

TABLE 3

Table 3. Hierarchical analysis of molecular variance (AMOVA).

Pairwise F_ST comparisons between populations showed that values ranging from a maximum differentiation of 0.0368 between Em17 and En10 to a minimum differentiation value of 0.0107 exhibited between Em25 and H01 (Table 4). Despite the low pairwise F_ST values found in our transect, all of them were significant after Bonferroni correction (Table 4). In addition, pairwise D_ST comparisons between populations indicated a similar pattern of genetic differentiation, with the higher D_ST value between Em25 and Em17 and between Em17 and En10 and the lower one between Em25 and H01 (Table 4).

TABLE 4

Table 4. Pairwise F_ST (above diagonal) and D_ST (below diagonal) values along the transect, estimated by means of Weir and Cockerham methods 1984 and as the harmonic mean across loci, respectively (Weir and Cockerham, 1984).

Genetic Structure and Gene Flow

The Bayesian inference of the genetic structure assigned the studied populations to two genetic clusters (K = 2), the second-most probable model being the one considering four genetic clusters (K = 4) (Supplementary Table 4). Considering K = 2, the plants belonging to the populations Em25, En10, and H01 exhibited very high membership proportions to a given ancestral genetic cluster, while the plants belonging to Em17 and En11 showed higher or medium values of membership proportions to the second cluster (Figure 3). However, when we consider the model assuming K = 4, the plants from E. mediohispanicum and E. nevadense populations showed medium to high membership proportions to four ancestral genetic clusters (Figure 3). In contrast, the individuals from the hybrid population exhibited medium to high values for their assignment probabilities to the most frequent genetic clusters in Em25 and En10 (Figure 3).

FIGURE 3

Figure 3. Significant gene flow values among E. mediohispanicum and E. nevadense populations surrounding the hybrid zone. Bottom: diagram with individual barplots from STRUCTURE (K = 2 and K = 4) based on the variation of 10 unlinked microsatellite loci.

The Bayesian inference indicated that the highest gene flow levels were taking place at the intrapopulational level for all the populations included in the study (Table 5). The gene flow among populations was not significant, except for the gene flow from Em25, En11, and En10 to the hybrid population (Table 5), the gene flow being higher from Em25 (29%) and En10 (28%) to H01 than that from En11 (14%).

TABLE 5

Table 5. Recent migration rates (m_ij) and SEs from i populations to j populations in the transect between E. mediohispanicum and E. nevadense, inferred from the variation at 10 microsatellite DNA loci using BayesAss.

Reproductive Barriers

Significant barriers to reproduction were found between E. mediohispanicum and E. nevadense, considering the population and species level. However, these barriers were asymmetric. The hybrid crosses performed on E. mediohispanicum plants presented a similar fitness to the intraspecific crosses performed in this species, considering both the cumulative predispersal component of the fitness and the total cumulative fitness (Table 6). When E. nevadense plants received the pollen of E. mediohispanicum, we found significant positive HI values at both the population and species levels, independent of the predispersal or cumulative fitness. Hence, hybrid crosses showed a mean decrease in fitness between 36% (considering Wpre in En11) and 56% (calculated for Wtot for En10). Considering the species level, the E. nevadense hybrids showed a mean decrease of 40% compared with the intraspecific crosses (Table 6).

TABLE 6

Table 6. Reproductive barrier estimation between E. mediohispanicum and E. nevadense.

Discussion

The consequences of new contacts between the closely related species have been widely highlighted for their interest in understanding the evolutionary processes (Barton and Hewitt, 1985; Harrison, 1993). More recently, they have been pointed out for their importance in assisting the adoption of better conservation policies in the face of rising hybridization (Gómez et al., 2015). In this sense, E. mediohispanicum and E. nevadense were suggested to be the two non-sister species (Abdelaziz et al., 2014) with contrasting distribution and a karyotypic history in the Iberian Peninsula (Muñoz-Pajares et al., 2017, 2018). However, as in the Sierra Nevada National Park, E. mediohispanicum inhabits the low lands, while E. nevadense has been described as an endemic species of the Sierra Nevada peaks (Blanca et al., 2009); both species could be the result of a recent or ongoing process of evolutionary divergence by the adaptation to high mountains. This elevational parapatric distribution could promote divergent ecological selection patterns, as those observed in other plant species (Campbell, 2004; Muir et al., 2013), preventing gene flow into the parental areas but creating a stable hybrid zone at the contact. This is supported in our case by the phenotypic stability of the contact zone after 10 years, and it is especially relevant because it makes the scenarios of non-sister species and non-completed speciation compatible as the origin of this hybrid zone. These hybridizing contact zones have important effects on the evolutionary dynamics of the hybridizing species and the biodiversity in the area where the contact happens (Wielstra, 2019; Alves de Moura et al., 2020). The hybridization is not always evident for human, which could have even more pervasive consequences for the biology of the species and their conservation (Gómez et al., 2015; Waldron et al., 2019).

With almost 200 genotyped plants along a transect from E. mediohispanicum to E. nevadense, we found significant differences between the populations, except for the putative hybrid population. This population seemed to be a mixture between genotypes coming from Em25 and En10, the populations in the extremes of the transect. Most of the genotypes of the hybrid population showed variable admixture proportions with the main genetic clusters of populations Em25 and En10. However, there were also plant genotypes belonging to one of the two inferred genetic sources in Em25 and En10 (Figure 3). Therefore, this hybrid zone shows the typical pattern found in the unimodal hybrid zones, with the probable presence of parental plants and/or backcrosses among hybrids, the more frequent plants in the population. These unimodal hybrid zones are frequent when the hybridizing species have not developed strong reproductive barriers (Harrison and Bogdanowicz, 1997). However, Jiggins and Mallet (2000) proposed that the unimodal and bimodal hybrid zones represent different stages of the speciation process: the former would correspond to the early stages in the speciation continuum, while the latter would signify that speciation is nearly completed (Mallet and Dasmahapatra, 2012).

The absence of an interspecific genetic structure and differentiation found in the altitudinal transect from E. mediohispanicum to E. nevadense seems to be related to the significant and unidirectional levels of gene flow detected for the hybrid population. These patterns are probably related to the highly diverse pollinator assemblage interacting with these species in the Sierra Nevada (Gómez et al., 2009a,b; Abdelaziz et al., unpublished data) and to their ability to promote pollen movement between populations (Tochigi et al., 2021; Abdelaziz et al., unpublished data), as supported by the high number of multilocus genotypes identified in this study. We could identify the hybrid zone despite the nearby distribution of our populations in this area. Usually, the hybrid zones have been studied in extended areas where the researchers identify parental populations and hybrids in a latitudinal gradient (e.g., Lepais et al., 2009; Liu et al., 2014). However, when searching for hybrid contact zones on altitudinal gradients, the approach must be more local: first, because the area at the mountain peaks is reduced with the elevation, and second, because the changes of kilometers in latitude correspond to the changes of meters in elevation for the environmental conditions affecting the plants (Jump et al., 2009). The interest for studying hybrid zones in elevation gradients is increasing in recent years (Aizawa and Iwaizumi, 2020; Mimura and Suga, 2020; Tamaki and Yamada, 2020), although some of these hybrid zones are already in the central body knowledge of hybridization in plants (Aguilar et al., 1999; James and Abbott, 2005).

The gene flow associated with the seed or pollen movements is not enough to generate hybrids in a contact zone. Weak reproductive barriers between the contacting species are also needed, whatever the mechanisms or the way to estimate those (Baack et al., 2015). We calculated pre- and post-dispersal components of reproductive barriers between E. mediohispanicum and E. nevadense. All of them were coherent and found evidence of the existence of asymmetric reproductive barriers between both species. In this sense, the hybrids obtained by pollinating the plants of E. nevadense experienced significant reductions of fitness when compared with the intraspecific crosses at both the predispersal (seed production) and postdispersal (germination and survival) components. In contrast, the hybrid crosses performed on E. mediohispanicum showed no fitness differences with the intraspecific crosses. Moreover, the E. mediohispanicum hybrids displayed a hybrid vigor for such traits as the number of stalks and the number of flowers, while significantly reduced values for these traits were found for E. nevadense hybrids. This means that the E. nevadense hybrids suffer from an important reduction in fitness at the prezygotic and postzygotic levels. The traits found to be modified in hybrids are related to floral display and are far from irrelevant. Indeed, they were recurrently demonstrated to play an important role in the attraction of pollinators and being under significant natural selection regimes in Erysimum (Gómez et al., 2009a) and other plant species (e.g., Johnston, 1991; Parachnowitsch and Kessler, 2010; Trunschke et al., 2017). Hence, the hybrids of E. mediohispanicum probably play an important role in maintaining the hybrid population.

In nature, the plants in the hybrid zone showed vigor only for traits related to the plant size (number of stalks and flowers), but intermediate values were shown for the rest (including all the flower size traits). Reductions in fitness have been described for hybrids mainly associated with sterility and developmental problems (Bomblies and Weigel, 2007; Ouyang et al., 2010; Levin, 2012), promoting the rise of efficient reproductive barriers (Ispolatov and Doebeli, 2009). However, hybrid organisms are also more vigorous than either of the parentals due to heterosis (Huang et al., 2015) or transgressive segregation in the further hybrid generations (Rieseberg et al., 1996, 1999). However, we have found that, while some traits can show heterotic patterns, some other traits can show just intermediate values to both parents. Depending on the trait modified by the hybridization, adaptive trajectories can significantly change and the evolutionary path of the hybrid population can move closer to or further away from the parental species. This fact, together with the absence of the important differences in flower size, may help in the occurrence of hybrids and backcrosses in the natural populations, favoring the continuity or even the spread of the hybrid population. Our phenotypic analysis, 10 years after the first sampling, confirms the same phenotypic pattern found in the transect in the first year and also found in the artificial hybrids developed in the greenhouse. Therefore, we can say that the hybrids are still produced in this contact zone. Studies on the stability of hybrids beyond F1 and on the natural selection acting on the parental and hybrid plants would clarify the stability of this hybrid zone and its effects on the biodiversity of the Sierra Nevada high mountains.

Despite the amount of data included in the current study, we could neither identify the reason why most of the gene flow to the hybrid population is coming from the most distant parental populations nor identify the causes of the asymmetry in the hybrid fitness. On the one hand, the non-homogeneous pattern of gene flow found in our altitudinal gradient could be related to the erratic movements or long-distance flights of the pollinators, as those described for Syzygium tierneyanum and Delphinium nuttallianum (Hopper, 1980; Schulke and Waser, 2001). The composition of the pollinator assemblage has an important impact on the gene flow and the genetic structure pattern found on the plant species interacting with them (Brunet et al., 2019). On the other hand, the mechanism underlying the asymmetry in hybrid fitness could be related to asymmetries in the interaction between the pollen and the stigma. A similar interaction seems to play a central role in the asymmetric patterns of hybridization found in different plant species (Rahmé et al., 2009; Pickup et al., 2019). Moreover, the outbreeding depression (Barmentlo et al., 2018) promoted by the intense patterns of local adaptation to the high elevation of the mountain could significantly contribute to the asymmetrical hybridization pattern found in this study.

Conclusions

The highly coherent results obtained by combining the field, the experimental, and the genetic data demonstrate that, in the altitudinal gradient between E. mediohispanicum and E. nevadense studied in the present work, there exists a narrow unimodal hybrid zone. The hybrid population appears to be maintained by high levels of the incoming gene flow from mainly two parental populations and by the weak and asymmetrical reproductive barriers between the parental species. E. mediohispanicum tolerates the interspecific pollen, but this does not occur with E. nevadense. This places E. mediohispanicum in a better situation for introgression events than E. nevadense. This pattern was replicated in the greenhouse, where E. mediohispanicum hybrids presented some heterotic traits, while E. nevadense showed reduced values for these same traits. Finally, a similar pattern was also found in the field when the same traits related to plant sizes showed heterotic values in the hybrid population compared with those in the parental populations. The heterotic phenotypes are probably more frequent in the hybrid zone because of their advantage in attracting pollinators and improving fitness since these traits were demonstrated to be under selection in E. mediohispanicum (Gómez et al., 2009a). Finding the same pattern after 10 years is evidence of the stability of this hybrid zone between E. mediohispanicum and E. nevadense at the top of the Sierra Nevada mountains.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

MA, AJM-P, FP, and JMG conceived and designed the experiments. MA, AJM-P, MB, and AG-M performed the experiments. MA, AJM-P, FP, JMG, MB, and AG-M analyzed the data. MA, AJM-P, FP, and JMG drafted and modified the manuscript. All authors read and approved the manuscript.

Funding

This research has been supported by a grant from the Spanish Ministry of Economy and Competitiveness (CGL2014-59886-JIN), the Organismo Autónomo de Parques Nacionales (Ref: 2415/2017), and the Ministry of Science and Innovation (PID2019-111294GB-I00/SRA/10.13039/501100011033). FP acknowledges the projects A.RNM.505.UGR18 (FEDER/Junta de Andalucía-Consejería de Economía y Conocimiento) and the Spanish Ministry of Science and Competitiveness (CGL2016-79950-R; CGL2017-86626-C2-2-P), including FEDER funds. MA was supported by the TransSpeciation project (CGL2014-59886-JIN). AJM-P was funded by the European Commission under the Marie Sklodowska-Curie Action Cofund 2016 EU agreement 754446 and the UGR Research and knowledge Transfer—Athenea3i. MB was supported by globalHybrids project (Ref: 2415/2017) and AG-M was supported by OUTevolution project (PID2019-111294GB-I00/SRA/10.13039/501100011033).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors thank the Sierra Nevada National Park headquarters for providing the permits to work into the National Park. The authors are grateful to Belén Herrador for her assistance during the lab work. The authors also thank the comments on their work of the associated editor and two reviewers, which certainly improved the previous version of this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.687094/full#supplementary-material

References

Ågren, J., and Schemske, D. W. (1993). Outcrossing rate and inbreeding depression in two annual monoecious herbs, Begonia hirsuta and B. semiovata. Evolution 47, 125–135. doi: 10.1111/j.1558-5646.1993.tb01204.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Abbott, R. J. (2017). Plant speciation across environmental gradients and the occurrence and nature of hybrid zones. J. Syst. Evol. 55, 238–258. doi: 10.1111/jse.12267

CrossRef Full Text | Google Scholar

Abdelaziz, M., Lorite, J., Muñoz-Pajares, A. J., Herrador, M. B., Perfectti, F., and Gómez, J. M. (2011). Using complementary techniques to distinguish cryptic species: a new erysimum (Brassicaceae) species from North Africa. Am. J. Bot. 98, 1049–1060. doi: 10.3732/ajb.1000438

PubMed Abstract | CrossRef Full Text | Google Scholar

Abdelaziz, M., Muñoz-Pajares, A. J., Lorite, J., Herrador, M. B., Perfectti, F., and Gómez, J. M. (2014). Phylogenetic relationships of Erysimum (Brassicaceae) from the baetic mountains (SE Iberian Peninsula. Anal. Jardín Botán. Madrid 1, 1–10. doi: 10.3989/ajbm.2377

CrossRef Full Text | Google Scholar

Aguilar, J. F., Rossell,ó, J. A, and Feliner, G. N. (1999). Molecular evidence for the compilospecies model of reticulate evolution in Armeria (Plumbaginaceae). Syst. Biol. 48, 735–754. doi: 10.1080/106351599259997

PubMed Abstract | CrossRef Full Text | Google Scholar

Aizawa, M., and Iwaizumi, M. G. (2020). Natural hybridization and introgression of Abies firma and Abies homolepis along the altitudinal gradient and genetic insights into the origin of Abies umbellata. Plant Species Biol. 35, 147–157. doi: 10.1111/1442-1984.12269

CrossRef Full Text | Google Scholar

Al-Shebaz, I. A. (2012). A generic and tribal synopsis of the Brassicaceae (Cruciferae). Taxon 6, 931–954. doi: 10.1002/tax.615002

CrossRef Full Text | Google Scholar

Alves de Moura, Y., Alves-Pereira, A., Silva, C. C., de Souza, L. M., de Souza, A. P., and Koehler, S. (2020). Secondary origin, hybridization and sexual reproduction in a diploid-tetraploid contact zone of the facultatively apomictic orchid Zygopetalum mackayi. Plant Biol. 22, 939–948. doi: 10.1111/plb.13148

PubMed Abstract | CrossRef Full Text | Google Scholar

Ancev, M. (2006). Polyploidy and hybridization in bulgarian Brassicaceae: distribution and evolutionary role. Phytol. Balcan. 12, 357–366.

Google Scholar

Anderson, E. (1949). Introgressive Hybridization. New York, NY: John Wiley. doi: 10.5962/bhl.title.4553

CrossRef Full Text | Google Scholar

Baack, E., Melo, M. C., Rieseberg, L. H., and Ortiz-Barrientos, D. (2015). The origins of reproductive isolation in plants. New Phytol. 207, 968–984. doi: 10.1111/nph.13424

PubMed Abstract | CrossRef Full Text | Google Scholar

Barmentlo, S. H., Meirmans, P. G., Luijten, S. H., Triest, L., and Oostermeijer, J. G. B. (2018). Outbreeding depression and breeding system evolution in small, remnant populations of Primula vulgaris: consequences for genetic rescue. Conserv. Genet. 19, 545–554. doi: 10.1007/s10592-017-1031-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Barton, N. H., and Hewitt, G. M. (1985). Analysis of hybrid zones. Annu. Rev. Ecol. Evol. Syst. 16, 113–148. doi: 10.1146/annurev.es.16.110185.000553

CrossRef Full Text | Google Scholar

Blanca, G., Cabezudo, B., Cueto, M., López, C. F., and Torres, C. M. (2009). Flora Vascular de Andalucía Oriental. Consejería de Medio Ambiente, Junta de Andalucía, Seville, Spain.

Google Scholar

Bomblies, K., and Weigel, D. (2007). Hybrid necrosis: autoimmunity as a potential gene-flow barrier in plant species. Nat. Rev. Genet. 8, 382–393. doi: 10.1038/nrg2082

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunet, J., Zhao, Y., and Clayton, M. K. (2019). Linking the foraging behavior of three bee species to pollen dispersal and gene flow. PLoS ONE 14:e0212561. doi: 10.1371/journal.pone.0212561

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, D., and Aldridge, G. (2006). “Floral biology of hybrid zones,” in Ecology and Evolution of Flowers, eds L. D. Harder, and S. C. H. Barrett (Oxford: Oxford University Press).

Google Scholar

Campbell, D. R. (2004). Natural selection in Ipomopsis hybrid zones: implications for ecological speciation. New Phytol. 161, 83–90. doi: 10.1046/j.1469-8137.2003.00919.x

CrossRef Full Text | Google Scholar

Canty, A., and Ripley, B. (2017). Package ‘Boot'. Bootstrap Functions. Version. 1.3.

Clot, B. (1992). Caryosystématique de quelques Erysimum L. dans le nord de la Péninsule Ibérique. Anal. Jardín Botán. Madrid 49, 215–229.

Google Scholar

Couvreur, T. L. P., Franzke, A., Al-Shehbaz, I. A., Bakker, F. T., Koch, M. A., and Mummenhoff, K. (2010). Molecular phylogenetics, temporal diversification, and principles of evolution in the mustard family (Brassicaceae). Mol. Biol. Evol. 27, 55–71. doi: 10.1093/molbev/msp202

PubMed Abstract | CrossRef Full Text | Google Scholar

Coyne, J. A., and Orr, H. A. (2004). Speciation. Sunderland, MA: Sinauer Ass., Inc.

Google Scholar

Cruz, R., Rolán-Álvarez, E., and García, C. (2001). Sexual selection on phenotypic traits in a hybrid zone of Littorina saxatilis (olivi). J. Evol. Biol. 14, 773–785. doi: 10.1046/j.1420-9101.2001.00324.x

CrossRef Full Text | Google Scholar

Cruzan, M. B., and Arnold, M. L. (1993). Ecological and genetic associations in an Irish hybrid zone. Evolution 47, 1432–1445. doi: 10.1111/j.1558-5646.1993.tb02165.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobzhansky, T. (1940). Speciation as a stage in evolutionary divergence. Am. Nat. 74, 312–321. doi: 10.1086/280899

CrossRef Full Text | Google Scholar

Dowling, T. E., and Secor, C. L. (1997). The role of hybridization and introgression in the diversification of animals. Annu. Rev. Ecol. Syst. 28, 593–619. doi: 10.1146/annurev.ecolsys.28.1.593

CrossRef Full Text | Google Scholar

Earl, D. A., and von Holdt, B. M. (2012). Structure harvester: a website and program for visualizing STRUCTURE output and implementing the evanno method. Conserv. Genet. Resour. 4, 359–361. doi: 10.1007/s12686-011-9548-7

CrossRef Full Text | Google Scholar

Ellstrand, N. C., Whitkus, R., and Rieseberg, L. H. (1996). Distribution of spontaneous plant hybrids. Proc. Natl. Acad. Sci. U.S.A. 93, 5090–5093. doi: 10.1073/pnas.93.10.5090

PubMed Abstract | CrossRef Full Text | Google Scholar

Emms, S. K., and Arnold, M. L. (1997). The effect of habitat on parental and hybrid fitness: transplant experiments with Louisiana irises. Evolution 51, 1112–1119. doi: 10.1111/j.1558-5646.1997.tb03958.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Evanno, G., Regnaut, S., and Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Excoffier, L., and Lischer, H. E. L. (2010). Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under linux and windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Falush, D., Stephens, M., and Pritchard, J. K. (2003). Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164, 1567–1587. doi: 10.1093/genetics/164.4.1567

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuertes Aguilar, J., Rossell,ó, J. A., and Nieto Feliner, G. (1999). Nuclear ribosomal DNA (nrDNA) concerted evolution in natural and artificial hybrids of Armeria (Plumbaginaceae). Mol. Ecol. 8, 1341–1346. doi: 10.1046/j.1365-294X.1999.00690.x

PubMed Abstract | CrossRef Full Text

Gómez, J. M. (2003). Herbivory reduces the strength of pollinator-mediated selection in the Mediterranean herb Erysimum mediohispanicum: consequences for plant specialization. Am. Nat. 162, 242–256. doi: 10.1086/376574

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez, J. M., Abdelaziz, M., Camacho, J. P. M., Muñoz-Pajares, A. J., and Perfectti, F. (2009b). Local adaptation and maladaptation to pollinators in a generalist geographic mosaic. Ecol. Lett. 12, 672–682. doi: 10.1111/j.1461-0248.2009.01324.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez, J. M., Bosch, J., Perfectti, F., Fernández, J., and Abdelaziz, M. (2007). Pollinator diversity affects plant reproduction and recruitment: the tradeoffs of generalization. Oecologia 153, 597–605. doi: 10.1007/s00442-007-0758-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez, J. M., González-Megías, A., Lorite, J., Abdelaziz, M., and Perfectti, F. (2015). The silent extinction: climate change and the potential hybridization-mediated extinction of endemic high-mountain plants. Biodivers. Conserv. 24, 1843–1857. doi: 10.1007/s10531-015-0909-5

CrossRef Full Text | Google Scholar

Gomez, J. M., Munoz-Pajares, A. J., Abdelaziz, M., Lorite, J., and Perfectti, F. (2014). Evolution of pollination niches and floral divergence in the generalist plant Erysimum mediohispanicum. Ann. Bot. 113, 237–249. doi: 10.1093/aob/mct186

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez, J. M., Perfectti, F., Bosch, J., and Camacho, J. P. M. (2009a). A geographic selection mosaic in a generalized plant–pollinator–herbivore system. Ecol. Monogr. 79, 245–264. doi: 10.1890/08-0511.1

CrossRef Full Text | Google Scholar

Good, T. P., Ellis, J. C., Annett, C. A., and Pierotti, R. (2000). Bounded hybrid superiority in an avian hybrid zone: effects of mate, diet, habitat choice. Evolution 54, 1774–1783. doi: 10.1111/j.0014-3820.2000.tb00721.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Goudet, J. (2005). HIERFSTAT, a package for R to compute and test hierarchical F-statistics. Mol. Ecol. Notes 5, 184–186. doi: 10.1111/j.1471-8286.2004.00828.x

CrossRef Full Text | Google Scholar

Greuter, W., Burdet, H. M., and Long, G. (1986). Med-checklist 3, Dicotyledones (Convolvulaceae-Labiatae. Gen è ve: Conservatoire et Jardin botaniques de la Ville de Gen è ve

Google Scholar

Harrison, R. G. (1993). Hybrid Zones and the Evolutionary Process. New York, NY: Oxford University Press.

Google Scholar

Harrison, R. G., and Bogdanowicz, S. M. (1997). Patterns of variation and linkage disequilibrium in a field cricket hybrid zone. Evolution 51, 493–505. doi: 10.1111/j.1558-5646.1997.tb02437.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hopper, S. D. (1980). Pollination of the rain-forest tree Syzygium tierneyanum (Myrtaceae) at Kuranda, Northern Queensland. Aust. J. Bot. 28, 223–237. doi: 10.1071/BT9800223

CrossRef Full Text | Google Scholar

Huang, X., Yang, S., Gong, J., Zhao, Y., Feng, Q., Gong, H., et al. (2015). Genomic analysis of hybrid rice varieties reveals numerous superior alleles that contribute to heterosis. Nat. Commun. 6:6258. doi: 10.1038/ncomms7258

PubMed Abstract | CrossRef Full Text

Ispolatov, I., and Doebeli, M. (2009). Speciation due to hybrid necrosis in plant–pathogen models. Evolut. Int. J. Organ. Evolut. 63, 3076–3084. doi: 10.1111/j.1558-5646.2009.00800.x

PubMed Abstract | CrossRef Full Text | Google Scholar

James, J. K., and Abbott, R. J. (2005). Recent, allopatric, homoploid hybrid speciation: the origin of Senecio squalidus (Asteraceae) in the British Isles from a hybrid zone on Mount Etna, Sicily. Evolution 59, 2533–2547. doi: 10.1111/j.0014-3820.2005.tb00967.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiggins, C. D., and Mallet, J. (2000). Bimodal hybrid zones and speciation. Trends Ecol. Evol. 15, 250–255. doi: 10.1016/S0169-5347(00)01873-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnston, M. O. (1991). Natural selection on floral traits in two species of Lobelia with different pollinators. Evolution 45, 1468–1479. doi: 10.1111/j.1558-5646.1991.tb02649.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jost, L. (2008). Gst and its relatives do not measure differentiation. Mol. Ecol. 17, 4015–4026. doi: 10.1111/j.1365-294X.2008.03887.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jump, A. S., Mátyás, C., and Peñuelas, J. (2009). The altitude-for-latitude disparity in the range retractions of woody species. Trends Ecol. Evolut. 24, 694–701. doi: 10.1016/j.tree.2009.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Keller, B., Ganz, R., Mora-Carrera, E., Nowak, M. D., Theodoridis, S., Koutroumpa, K., et al. (2021). Asymmetries of reproductive isolation are reflected in directionalities of hybridization: integrative evidence on the complexity of species boundaries. New Phytol. 229, 1795–1809. doi: 10.1111/nph.16849

PubMed Abstract | CrossRef Full Text | Google Scholar

Key, K. H. L. (1968). The concept of stasipatric speciation. Syst. Zool. 17, 14–22. doi: 10.1093/sysbio/17.1.14

CrossRef Full Text | Google Scholar

Knief, U., Bossu, C. M., Saino, N., Hansson, B., Poelstra, J., Vijay, N., et al. (2019). Epistatic mutations under divergent selection govern phenotypic variation in the crow hybrid zone. Nat. Ecol. Evolut. 3, 570–576. doi: 10.1038/s41559-019-0847-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Lepais, O., Petit, R. J., Guichoux, E., Lavabre, J. E., Alberto, F., Kremer, A., et al. (2009). Species relative abundance and direction of introgression in oaks. Mol. Ecol. 18, 2228–2242. doi: 10.1111/j.1365-294X.2009.04137.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Levin, D. A. (1967). Hybridization between anual species of phlox: population structure. Am. J. Bot. 54, 1122–1130. doi: 10.1002/j.1537-2197.1967.tb10742.x

CrossRef Full Text | Google Scholar

Levin, D. A. (2012). The long wait for hybrid sterility in flowering plants. New Phytol. 196, 666–670. doi: 10.1111/j.1469-8137.2012.04309.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, B., Abbott, R. J., Lu, Z., Tian, B., and Liu, J. (2014). Diploid hybrid origin of ostryopsis intermedia (Betulaceae) in the Qinghai-Tibet Plateau triggered by quaternary climate change. Mol. Ecol. 23, 3013–3027. doi: 10.1111/mec.12783

PubMed Abstract | CrossRef Full Text | Google Scholar

Mallet, J., and Dasmahapatra, K. K. (2012). Hybrid zones and the speciation continuum in Heliconius butterflies. Mol. Ecol. 21, 5643–5645. doi: 10.1111/mec.12058

PubMed Abstract | CrossRef Full Text | Google Scholar

Marhold, K., and Lihová, J. (2006). Polyploidy, hybridization and reticulate evolution: lessons from the Brassicaceae. Plant Syst. Evolut. 259, 143–174. doi: 10.1007/s00606-006-0417-x

CrossRef Full Text | Google Scholar

Marques, I., Feliner, G. N., Draper Munt, D., Martins-Loução, M. A., and Aguilar, J. F. (2010). Unraveling cryptic reticulate relationships and the origin of orphan hybrid disjunct populations in Narcissus. Evolut. Int. J. Organ. Evolut. 64, 2353–2368. doi: 10.1111/j.1558-5646.2010.00983.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayr, E. (1942). Systematics and the Origin of Species. New York, NY: Columbia University Press.

Google Scholar

McMillan, W. O., Jiggins, C. D., and Mallet, J. (1997). What initiates speciation in passion-vine butterflies?. Proc. Natl. Acad. Sci. U.S.A. 94, 8628–8633. doi: 10.1073/pnas.94.16.8628

PubMed Abstract | CrossRef Full Text | Google Scholar

Mimura, M., and Suga, M. (2020). Ambiguous species boundaries: Hybridization and morphological variation in two closely related rubus species along altitudinal gradients. Ecol. Evol. 10, 7476–7486. doi: 10.1002/ece3.6473

PubMed Abstract | CrossRef Full Text | Google Scholar

Minder, A. M., Rothenbuehler, C., and Widmer, A. (2007). Genetic structure of hybrid zones between Silene latifolia and Silene dioica (Caryophyllaceae): evidence for introgressive hybridization. Mol. Ecol. 16, 2504–2516. doi: 10.1111/j.1365-294X.2007.03292.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohamed, M. A. (2013). How species are evolutionarily maintained: Pollinator-mediated divergence and hybridization in Erysimum mediohispanicum and Erysimum nevadense (Doctoral dissertation). Universidad de Granada.

Google Scholar

Muir, G., Osborne, O. G., Sarasa, J., Hiscock, S. J., and Filatov, D. A. (2013). Recent ecological selection on regulatory divergence is shaping clinal variation in Senecio on Mount Etna. Evolution 67, 3032–3042. doi: 10.1111/evo.12157

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz-Pajares, A. J., Abdelaziz, M., and Picó, F. X. (2020). Temporal migration rates affect the genetic structure of populations in the biennial Erysimum mediohispanicum with reproductive asynchrony. AoB Plants 12:plaa037. doi: 10.1093/aobpla/plaa037

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz-Pajares, A. J., García, C., Abdelaziz, M., Bosch, J., Perfectti, F., and Gómez, J. M. (2017). Drivers of genetic differentiation in a generalist insect-pollinated herb across spatial scales. Mol. Ecol. 26, 1576–1585. doi: 10.1111/mec.13971

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz-Pajares, A. J., Herrador, M. B., Abdelaziz, M., Picó, F. X., Sharbel, T. F., Gómez, J. M., et al. (2011). Characterization of microsatellite loci in Erysimum mediohispanicum (Brassicaceae) and cross-amplification in related species. Am. J. Bot. 98, e287–e289. doi: 10.3732/ajb.1100181

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz-Pajares, A. J., Perfectti, F., Loureiro, J., Abdelaziz, M., Biella, P., Castro, M., et al. (2018). Niche differences may explain the geographic distribution of cytotypes in Erysimum mediohispanicum. Plant Biol. 20, 139–147. doi: 10.1111/plb.12605

PubMed Abstract | CrossRef Full Text | Google Scholar

Nei, M. (1987). Molecular Evolutionary Genetics. New York, NY: Columbia University Press. doi: 10.7312/nei-92038

CrossRef Full Text | Google Scholar

Nieto Feliner, G., Rosato, M., Alegre, G., San Segundo, P., Rossell,ó, J. A, Garnatje, T., et al. (2019). Dissimilar molecular and morphological patterns in an introgressed peripheral population of a sand dune species (Armeria pungens, Plumbaginaceae. Plant Biol. 21, 1072–1082. doi: 10.1111/plb.13035

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, Y., Liu, Y. G., and Zhang, Q. (2010). Hybrid sterility in plant: stories from rice. Curr. Opin. Plant Biol. 13, 186–192. doi: 10.1016/j.pbi.2010.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Owens, G. L., Baute, G. J., and Rieseberg, L. H. (2016). Revisiting a classic case of introgression: hybridization and gene flow in Californian sunflowers. Mol. Ecol. 25, 2630–2643. doi: 10.1111/mec.13569

PubMed Abstract | CrossRef Full Text | Google Scholar

Parachnowitsch, A. L., and Kessler, A. (2010). Pollinators exert natural selection on flower size and floral display in penstemon digitalis. New Phytol. 188, 393–402. doi: 10.1111/j.1469-8137.2010.03410.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Parsons, K. J., Son, Y. H., and Albertson, R. C. (2011). Hybridization promotes evolvability in African cichlids: connections between transgressive segregation and phenotypic integration. Evol. Biol. 38, 306–315. doi: 10.1007/s11692-011-9126-7

CrossRef Full Text | Google Scholar

Pickup, M., Brandvain, Y., Fraïsse, C., Yakimowski, S., Barton, N. H., Dixit, T., et al. (2019). Mating system variation in hybrid zones: facilitation, barriers and asymmetries to gene flow. New Phytol. 224, 1035–1047. doi: 10.1111/nph.16180

PubMed Abstract | CrossRef Full Text | Google Scholar

Polatschek, A. (1986). “Erysimum,” in Mountain flora of Greece, Vol. l, ed A. Strid (Cambridge: Cambridge University Press), 239–247.

Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959. doi: 10.1093/genetics/155.2.945

PubMed Abstract | CrossRef Full Text | Google Scholar

R Development Core Team. (2014). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Available online at: http://www.R-project.org

Google Scholar

Rahmé, J., Widmer, A., and Karrenberg, S. (2009). Pollen competition as an asymmetric reproductive barrier between two closely related Silene species. J. Evol. Biol. 22, 1937–1943. doi: 10.1111/j.1420-9101.2009.01803.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rambaut, A., and Drummond, A. J. (2007). Tracer v1.4 [Computer Program]. Available online at: http://tree.bio.ed.ac.uk/software/tracer/

Rieseberg, L. H., Archer, M. A., and Wayne, R. K. (1999). Transgressive segregation, adaptation and speciation. Heredity 83, 363–372. doi: 10.1038/sj.hdy.6886170

PubMed Abstract | CrossRef Full Text | Google Scholar

Rieseberg, L. H., Baird, S. J. E., and Desroche, A. (1998). Patterns of mating in wild sunflower hybrid zones. Evolution 52, 713–726. doi: 10.1111/j.1558-5646.1998.tb03696.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rieseberg, L. H., Sinervo, B., Linder, C. R., Ungerer, M. C., and Arias, D. M. (1996). Role of gene interactions in hybrid speciation: evidence from ancient and experimental hybrids. Science 272, 741–745. doi: 10.1126/science.272.5262.741

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulke, B., and Waser, N. M. (2001). Long-distance pollinator flights and pollen dispersal between populations of Delphinium nuttallianum. Oecologia 127, 239–245. doi: 10.1007/s004420000586

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamaki, I., and Yamada, Y. (2020). Environmental pressure rather than ongoing hybridization is responsible for an altitudinal cline in the morphologies of two oaks. J. Plant Ecol. 13, 413–422. doi: 10.1093/jpe/rtaa028

CrossRef Full Text | Google Scholar

Tastard, E., Ferdy, J. B., Burrus, M., Thébaud, C., and Andalo, C. (2012). Patterns of floral colour neighbourhood and their effects on female reproductive success in an antirrhinum hybrid zone. J. Evol. Biol. 25, 388–399. doi: 10.1111/j.1420-9101.2011.02433.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tochigi, K., Shuri, K., Kikuchi, S., Naoe, S., Koike, S., and Nagamitsu, T. (2021). Phenological shift along an elevational gradient and dispersal of pollen and seeds maintain a hybrid zone between two cherry tree species. Plant Species Biol. 36, 230–245. doi: 10.1111/1442-1984.12311

CrossRef Full Text | Google Scholar

Trunschke, J., Sletvold, N., and Ågren, J. (2017). Interaction intensity and pollinator-mediated selection. New Phytol. 214, 1381–1389. doi: 10.1111/nph.14479

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, B. L. (2006). Taxonomy and nomenclature of the Erysimum asperum- E. capitatum complex (Brassicaceae). Phytologia 88, 279–287. doi: 10.5962/bhl.part.10454

CrossRef Full Text | Google Scholar

Vedenina, V. Y., and Helversen, O. V. (2003). Complex courtship in a bimodal grasshopper hybrid zone. Behav. Ecol. Sociobiol. 54, 44–54. doi: 10.1007/s00265-003-0595-2

CrossRef Full Text | Google Scholar

Waldron, B. P., Kuchta, S. R., Hantak, M. M., Hickerson, C. A. M., and Anthony, C. D. (2019). Genetic analysis of a cryptic contact zone between mitochondrial clades of the eastern red-backed salamander, Plethodon cinereus. J. Herpetol. 53, 144–153. doi: 10.1670/18-088

CrossRef Full Text | Google Scholar

Weir, B. S., and Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution 38, 1358–1370. doi: 10.1111/j.1558-5646.1984.tb05657.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitham, T. G., Martinsen, G. D., Keim, P., Floate, K. D., Dungey, H. S., and Potts, B. M. (1999). Plant hybrid zones affect biodiversity: tools for a genetic-based understanding of community structure. Ecology 80, 416–428. doi: 10.1890/0012-9658(1999)0800416:PHZABT2.0.CO;2

CrossRef Full Text | Google Scholar

Whitham, T. G., Morrow, P. A., and Potts, B. M. (1994). Plant hybrid zones as centers of biodiversity: the herbivore community of two endemic tasmanian eucalypts. Oecologia 97, 481–490. doi: 10.1007/BF00325886

PubMed Abstract | CrossRef Full Text | Google Scholar

Wielstra, B. (2019). Historical hybrid zone movement: More pervasive than appreciated. J. Biogeogr. 46, 1300–1305. doi: 10.1111/jbi.13600

CrossRef Full Text | Google Scholar

Wielstra, B., Burke, T., Butlin, R. K., Avci, A., Üzüm, N., Bozkurt, E, et al. (2017). A genomic footprint of hybrid zone movement in crested newts. Evolut. Lett. 1, 93–101. doi: 10.1002/evl3.9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, G. A., and Rannala, B. (2003). Bayesian inference of recent migration rates using multilocus genotypes. Genetics 163, 1177–1191. doi: 10.1093/genetics/163.3.1177

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, N. D. (1996). An analysis of the causes of genetic isolation in two Pacific coast iris hybrid zones. Can. J. Bot. 74, 2006–2013. doi: 10.1139/b96-241

CrossRef Full Text | Google Scholar