- 1Departamento de Genética, Facultad de Ciencias, Campus Fuentenueva, Universidad de Granada, Granada, Spain
- 2Laboratório Associado, Plant Biology, Research Centre in Biodiversity and Genetic Resources, Centro de Investigação em Biodiversidade e Recursos Genéticos, Universidade Do Porto, Campus Agrário de Vairão, Fornelo e Vairão, Portugal
- 3Research Unit Modeling Nature, Universidad de Granada, Granada, Spain
- 4Departamento de Ecología Funcional y Evolutiva, Estación Experimental de Zonas Áridas, Consejo Superior de Investigaciones Científicas, Almeria, Spain
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. (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. 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 (NG), as the number of genotypes showing at least one different allele, excluding missing data; (b) mean number of alleles per locus (na); (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 (RP), 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 (FIS), which provides information about the Hardy–Weinberg equilibrium departures due to either excess or defect of heterozygotes. We estimated FIS 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. 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 DST (Jost, 2008). FST 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 (WPRE), as fruit set x seed production, and the cumulative total fitness (WTOT), 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
where wo is the fitness of the intraspecific outcrossing treatments and wh 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. 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 (na) 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 HO 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 HS presented the same pattern shown by HO but the values of HS were higher than Ho in all the considered populations and species (Table 2). The values of the mean allelic richness (RS) 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 (RP) 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 (FIS) 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).
Pairwise FST 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 FST values found in our transect, all of them were significant after Bonferroni correction (Table 4). In addition, pairwise DST comparisons between populations indicated a similar pattern of genetic differentiation, with the higher DST value between Em25 and Em17 and between Em17 and En10 and the lower one between Em25 and H01 (Table 4).
Table 4. Pairwise FST (above diagonal) and DST (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. 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. Recent migration rates (mij) 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).
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
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
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
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
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
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
Al-Shebaz, I. A. (2012). A generic and tribal synopsis of the Brassicaceae (Cruciferae). Taxon 6, 931–954. doi: 10.1002/tax.615002
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
Ancev, M. (2006). Polyploidy and hybridization in bulgarian Brassicaceae: distribution and evolutionary role. Phytol. Balcan. 12, 357–366.
Anderson, E. (1949). Introgressive Hybridization. New York, NY: John Wiley. doi: 10.5962/bhl.title.4553
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
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
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
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.
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
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
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).
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
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.
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
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
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
Dobzhansky, T. (1940). Speciation as a stage in evolutionary divergence. Am. Nat. 74, 312–321. doi: 10.1086/280899
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Harrison, R. G. (1993). Hybrid Zones and the Evolutionary Process. New York, NY: Oxford University Press.
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
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
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
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
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
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
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
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
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
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
Key, K. H. L. (1968). The concept of stasipatric speciation. Syst. Zool. 17, 14–22. doi: 10.1093/sysbio/17.1.14
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
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
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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
Nei, M. (1987). Molecular Evolutionary Genetics. New York, NY: Columbia University Press. doi: 10.7312/nei-92038
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Trunschke, J., Sletvold, N., and Ågren, J. (2017). Interaction intensity and pollinator-mediated selection. New Phytol. 214, 1381–1389. doi: 10.1111/nph.14479
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
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
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
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
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
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
Wielstra, B. (2019). Historical hybrid zone movement: More pervasive than appreciated. J. Biogeogr. 46, 1300–1305. doi: 10.1111/jbi.13600
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
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
Keywords: hybridization, Erysimum mediohispanicum, Erysimum nevadense, Sierra Nevada, phenotype, reproductive isolation
Citation: Abdelaziz M, Muñoz-Pajares AJ, Berbel M, García-Muñoz A, Gómez JM and Perfectti F (2021) Asymmetric Reproductive Barriers and Gene Flow Promote the Rise of a Stable Hybrid Zone in the Mediterranean High Mountain. Front. Plant Sci. 12:687094. doi: 10.3389/fpls.2021.687094
Received: 28 March 2021; Accepted: 15 July 2021;
Published: 25 August 2021.
Edited by:
Andrew A. Crowl, Duke University, United StatesReviewed by:
Adrian Christopher Brennan, Durham University, United KingdomSantiago Martín-Bravo, Universidad Pablo de Olavide, Spain
Copyright © 2021 Abdelaziz, Muñoz-Pajares, Berbel, García-Muñoz, Gómez and Perfectti. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mohamed Abdelaziz, bWFiZGVsYXppem0mI3gwMDA0MDt1Z3IuZXM=
†These authors have contributed equally to this work and share first authorship