Long-term persisting hybrid swarm and geographic difference in hybridization pattern: genetic consequences of secondary contact between two Vincetoxicum species (Apocynaceae–Asclepiadoideae)
© Li et al. 2016
Received: 23 September 2015
Accepted: 12 January 2016
Published: 22 January 2016
During glacial periods, glacial advances caused temperate plant extirpation or retreat into localized warmer areas, and subsequent postglacial glacial retreats resulted in range expansions, which facilitated secondary contact of previously allopatric isolated lineages. The evolutionary outcomes of secondary contact, including hybrid zones, dynamic hybrid swarm, and resultant hybrid speciation, depends on the strengths of reproductive barriers that have arisen through epistatic and pleiotropic effects during allopatric isolation. The aim of this study was to demonstrate refugia isolation and subsequent secondary contact between two perennial Asclepioid species and to assess the genetic consequences of the secondary contact. We modeled the range shift of two ecologically distinct Vincetoxicum species using the species distribution model (SDM) and assessed the genetic consequences of secondary contact by combining morphological and genetic approaches. We performed morphometric analysis (592 individuals) and examined 10 nuclear microsatellites (671 individuals) in V. atratum, V. japonicum, and putative hybrid populations.
Multivariate analysis, model-based Bayesian analysis, and non-model-based discriminant analysis of principal components confirmed the hybridization between V. atratum and V. japonicum. High pollen fertility and a lack of linkage disequilibrium suggested that the hybrid populations may be self-sustaining and have persisted since V. atratum and V. japonicum came into contact during the post-glacial period. Moreover, our findings show that the pattern of hybridization between V. atratum and V. japonicum is unidirectional and differs among populations. Geographically-isolated hybrid populations exist as genetically distinct hybrid swarms that consist of V. atratum-like genotypes, V. japonicum-like genotypes, or admixed genotypes. In addition, Bayesian-based clustering analysis and coalescent-based estimates of long-term gene flow showed patterns of introgressive hybridization in three morphologically ‘pure’ V. japonicum populations.
In this study, we demonstrated that climatic oscillations during the Quaternary period likely led to species range shift and subsequently secondary contact. Hybrid populations may be self-sustaining and have persisted since V. atratum and V. japonicum came into contact during the post-glacial period. Pattern of hybridization between V. atratum and V. japonicum is unidirectional and differs among populations. We concluded that these differences in the genetic consequences of secondary contact are caused by historical colonization processes and/or natural selection.
In response to climatic oscillations over the last 2 million years (Ma), the geographical distributions of temperate plants have been subject to multiple contractions and expansions . During glacial periods, glacial advances caused temperate plant extirpation or retreat into localized warmer areas, and subsequent postglacial glacial retreats resulted in range expansions, which facilitated secondary contact of previously allopatric isolated lineages [2–4]. The evolutionary outcomes of secondary contact, including hybrid zones, dynamic hybrid swarm, and resultant hybrid speciation , depend on the strengths of reproductive barriers .
Assuming that the reproductive isolation between two species coming into contact was incomplete, first generation hybrids may mate with parental species. Frequently backcrossing between hybrids and their parental species will lead introgression of alleles across species boundaries, thereby affecting species history and delimitation . The dynamics of introgression for neutral alleles are thought to be depends on the reproductive system of the interacting taxa [8–11], but also on relative population abundance . Genes are assumed to flow from large populations of one species into much smaller populations. Recent theoretical models predicted that the direction of introgression of neutral alleles between lineages should proceed from the locally established species towards the expanding congener [13, 14]. In addition, genotype combinations might be filtered by selection factors that vary with the surroundings (genotype-by-environment interactions) [15–22].
In recent years, an increasing number of studies have employed species distribution modeling (SDM) or ecological niche modeling (ENM) to support refugia isolation and subsequent range expansion and corroborate phylogenetic approaches. The aim of this study was to demonstrate refugia isolation and subsequent secondary contact by modeling range shifts. This study presents one of the first applications of the species distribution model (SDM) to explain the direction of historical gene flow and its effects on the genetic consequences of secondary contact by predicting potential glacial distribution and postglacial colonization processes.
The Japanese Archipelago is thought to have been free of massive ice sheets during the last glacial maximum (LGM) with the exception of some high altitude areas . Previous phylogenetic studies showed that the postglacial range expansion process may have caused secondary contact between formerly isolated lineages of Fagus crenata , Aesculus turbinata , Padus grayana, Euonymus oxyphyllus, Magnolia obovata, Carpinus tschonoskii, C. japonica, C. laxiflora , Platycrater arguta , Rubus palmatus, and R. grayanus . Two ecologically distinct species of Vincetoxicum (Apocynaceae-Asclepiadoideae), V. atratum (Bunge) Morren et Decaisne and V. japonicum Morren et Decne., provide a suitable study system. Because the different habitat preferences of the species may result in different responses to climate change [3, 29–31], V. atratum and V. japonicum may have retreated into separate refugia during the glacial period and subsequent postglacial expansion may have facilitated secondary contact, resulting in hybridization.
Vincetoxicum atratum and V. japonicum are perennial herbs. Vincetoxicum atratum is distributed across the Sino-Japanese region including the Japanese Archipelago, Korean Peninsula, and China, whereas V. japonicum is distributed in the Japanese Archipelago and the Korean Peninsula . In Japan, V. atratum is found in grasslands or meadows across northern to southern regions of the archipelago, while V. japonicum is found in sunny meadows or on rocky beaches near seashores of central and southern regions of the archipelago . Because the habitat area has diminished due to habitat destruction and degradation, the number and sizes of V. atratum populations have decreased rapidly. Morphologically, the two species can be distinguished based on specific characteristics. Vincetoxicum atratum has a purple corolla, which is pubescent outside and glabrous inside, while V. japonicum has a yellowish white corolla, which is glabrous outside and pubescent inside. Leaves of V. atratum are villous on both sides while those of V. japonicum are densely pubescent . The flowering periods of these species partially overlap in early summer. Some dipteran species are pollinators of both V. atratum and V. japonicum . The low level of genetic diversification in the chloroplast DNA (cpDNA) of the two species in contrast with the considerable morphological divergence suggests that the species have undergone rapid morphological divergence .
In this study, we examined past hybridization and introgression between V. atratum and V. japonicum by combining morphological and genetic approaches with paleoclimatic analysis. First, we examined the hybrid status of individuals from seven morphologically putative hybrid populations using a morphometric approach. Second, we examined the genetic structure of morphologically ‘pure’ V. atratum, V. japonicum, and putative hybrid populations. In particular, we examined deviation from linkage equilibrium for nuclear microsatellite markers because linkage disequilibrium would be expected in a hybrid zone as a result of ongoing gene flow between the parent species [35–37]. Third, we predicted species paleo-distributions in response to the Quaternary climate change to model range shifts of V. atratum and V. japonicum and to infer past opportunities for secondary contact between the species. In this report, we address the following three questions. (1) What are the genetic and phenotypic patterns of morphologically ‘pure’ V. atratum, V. japonicum, and putative hybrid populations? (2) What is the outcome of secondary contact between V. atratum and V. japonicum? (3) What is the genetic consequence of hybridization between the two species?
Study site and sampling
DNA extraction and nuclear microsatellite scoring
Total DNA was extracted from approximately 100 mg of leaf tissue based on the method of Maki, Horie and Yokoyama . Ten nuclear simple sequence repeat (nrSSR) loci including vinc5, vinc107, vinc118 , Vpy 012, Vpy 013, Vpy 016, Vpy 018, Vpy 022 , Vkat3, and Vkat2 (Yamashiro et al., unpublished,) were scored. These 10 loci were amplified in three multiplex reactions (Additional file 1: Table S2). Multiplex PCRs were performed in 3-μL volumes containing approximately 50 ng genomic DNA, 1 × Type-it Multiplex PCR Master Mix (QIAGEN), and 0.2 μM each of forward and reverse primers. The reaction parameters were an initial denaturation step at 95 °C for 5 min followed by 26 cycles at 95 °C for 30 s, 60 °C for 90 s, and 72 °C for 60 s, and a final step at 60 °C for 30 min. PCR products were run with an internal size standard GeneScan LIZ-600 (Applied Biosystems) on an ABI 3100 Genetic Analyzer (Applied Biosystems). Allele sizes were scored using GENEMAPPER v3.7 software (Applied Biosystems).
Description of the morphological characteristics scored for morphologically ‘pure’ V. atratum and V. japonicum plants and morphologically intermediate plants, and the results of one way ANOVA or Kruskal-Wallis ANOVA of morphologically ‘pure’ V. atratum and V. japonicum plants
Species F (dF = 1)
Maximum leaf blade length (cm)
Maximum leaf blade width (cm)
Petiole length (cm)
Stem length (cm)a
Hair density on adaxial leaf surface (cm)b
Hair density on abaxial leaf surface (cm)c
Number of flowers
Pedicel length (cm)
Pedicel width (cm)
Corolla diameter (cm)
Hair density on corolla outside (cm)d
Hair density on corolla inside (cm)e
Corona diameter (cm)
Gynostegium length (cm)
Gynostegium diameter (cm)
Corpusculum length (cm)
Pollinia length (cm)
For 13 continuous variables, significant differences between allopatric V. atratum and V. japonicum were tested using one-way analysis of variance (ANOVA). Pearson correlation coefficients were calculated to eliminate highly correlated (Pearson’s correlation > 0.9) characteristics from further analyses. A Kolmogorov-Smirnov test was used to determine whether the data met the assumption of normal distribution, while a Levene’s test was used to test for homoscedasticity. For four ordinal variables, Kruskal-Wallis one-way ANOVA was performed to estimate differences between V. atratum and V. japonicum. Principal coordinate analysis (PCO) was used to visualize morphological variation among morphologically ‘pure’ V. atratum, V. japonicum populations and morphologically hybrid populations. All statistical analyses were performed using R programming language .
Pollen viability estimation
Pollen fertility was examined in three populations of V. atratum (A01–A03) and V. japonicum (J01–J03) and in all putative hybrid populations following the procedure of .
Population genetic analysis
For each population, the mean number of alleles per locus (N a), observed heterozygosity (H O), and expected heterozygosity (H E) were calculated using the program GenAlEx version 6 . Allelic richness (AR) and private allele richness (PAR) were estimated using HP-RARE 1.1 software . Genetic differentiation in all pair wise combinations of the populations was estimated from the fixation index (F ST,  for nuclear simple sequence repeats (nrSSRs) with 1000 permutations using Arlequin version 3.5 software . Tests for deviation from the Hardy-Weinberg equilibrium (HWE) at each locus were conducted with Arlequin version 3.5 software using 10,000 permutations . Tests for the linkage disequilibrium (LD) of all combinations of the loci for each population were conducted using GENEPOP 4.2 software . LD can be induced by genetic admixture between divergent gene pools [35, 36] and will decay if gene flow occurs between the divergent gene pools , thus, if no pair wise LD was detected among loci, we would expect that the hybrid population would not be maintained by ongoing gene flow.
Bayesian admixture analysis
A Bayesian model-based approach implemented in the program STRUCTURE version 2.2.3  and TESS version 2.3 software  was used to assess the genetic delimitation between V. atratum and V. japonicum and to estimate the proportion of each individual genome originating from each of the parental populations. We assigned individuals having an assignment value (q) greater than 0.9 into one of the parental lineages and classified those with values less than 0.9 as hybrids. For analysis using STRUCTURE, different probable numbers of genetic clusters (K) were estimated for all samples under the admixture model with correlated allele frequencies. Ten replicate runs were performed for each value of K ranging from 1 to 10. Each run included a burn-in of 10,000 steps followed by 100,000 Markov chain Monte Carlo (MCMC) simulations. The optimal value of K was calculated using the STRUCTURE HARVESTER program . The results from the 10 replicates of the optimal value of K were averaged using the CLUMPP program  and were visualized using the DISTRUCT 1.1 program .
TESS version 2.3 was used to estimate population structure by incorporating geographical coordinates of individuals as prior information to detect discontinuities in allele frequencies . The analyses were performed using a conditional autoregressive (CAR) model with a burn-in of 10,000 iterations followed by 60,000 iterations for each K value (the number of clusters) between 2 and 10. The optimal value of K was determined based on the lowest value of the deviance information criterion (DIC). The results from the 10 replicates with the lowest DIC values were averaged using CLUMPP  and were visualized using DISTRUCT 1.1 as in the STRUCTURE analyses .
Discriminant analysis of principal components
To detect complex genetic structure, a non-model-based discriminant analysis of principal components (DAPC) using the adegenet version 1.2.8 package  in R  was conducted based on the SSR data set. The prior clusters were defined by the groups obtained from each population. DAPC first transforms genotype into uncorrelated components using principle component analysis (PCA) and then performs a discriminant analysis on the retained principle components (PCs).
Estimation of long-term gene flow
We divided studied populations based on the genetic constitution of each population defined by STRUCTURE, namely, ‘pure’ V. atratum (A01–A07 and A09), ‘pure’ V. japonicum (J01–J08), south group of putative hybrid population (H01–H05), north group of putative hybrid population (H06–H07) and introgressed V. japonicum (J09–J11). We estimated bidirectional long-term migration rates among the five categories based on 10 nrSSR loci with a coalescent-based model implemented in the LAMARC 2.1 program . We randomly picked genotype data for 20 individuals in each categories, because migrate history in additional individuals is probably similar to that of its cohorts. We estimated impacts of migration relative to mutation rate (M = m/μ) and the mutation-scaled effective population size (Θ = 4N e μ). N e is the effective population size, m is the migration rate between two populations, and μ is the mutation rate per generation at the locus considered. We used logarithmic prior distributions for migration rate (0.01, 3000), and effective population size (0.00001, 10). We conducted Bayesian analyses under the Brownian motion model . Search parameters consisted of 10 initial chains with 5,000 samples, a sampling interval of 30 (150,000 steps), and a burn-in of 10,000 samples; two final chains with 50,000 samples and a sampling interval of 30 (1,500,000 steps); and a burn-in of 200,000 samples. For each pair wise comparison (i → j), the number of immigrants from population i into population j (N e m (i → j)) was estimated from N e m = Θ i *M (i → j) ⁄ 4, where Θ represents the mutation scaled effective population size (Θ = 4N e μ). Results are reported as most probable estimates (MPE) with 95 % confidence intervals (CIs). We used Tracer v1.6  to confirm sampling adequacy and convergence of parameters. All parameter estimates were well supported, with ESS (effective sample size) values exceeding 200.
Species distribution modeling
We used species distribution modeling (SDM) to identify the range of suitable habitats of V. atratum and V. japonicum during three periods, namely, the present, the mid-Holocene (6,000 ybp), and the last glacial maximum (LGM; 21,000 ybp). The maximum entropy algorithm (MAXENT, Phillips and Dudík; 2008), based on presence-only data, and the associated environmental covariates were used to model species distribution. Species occurrence data were collected from the S-Net data portal (http://science-net.kahaku.go.jp/) and our sampling locations (Additional file 1: Table S3). We restricted the spatial extent of the modeled area to only Japanese Archipelago for both species. SDM was constructed using information for 125 V. atratum and 128 V. japonicum sites to predict the geographic distributions of each species (Additional file 1: Table S3; Additional file 1: Figure S1). Collection coordinates were plotted in Google Earth (http://earth.google.com) to confirm that each set represented a reasonable location. A total of 19 environmental variables (Additional file 1: Table S4) were downloaded from the WorldClim website (www.worldclim.org; . We extracted the values of each of the 19 environmental variables for all occurrence recorder points using the DIVA-GIS program  and examined pair wise correlations using R . After evaluation, variables with low pair wise Pearson correlation coefficients (r < 0.70) were used for subsequent analyses (Additional file 1: Table S3). We used the simulation summaries of the Community Climate System Model version 4 (CCSM4) to create mid-Holecene and LGM climate data projections. The resolutions of current, mid-Holecene, and LGM climate layers were 30 arcsec, 30 arcsec, and 2.5 arcmin, respectively. The distributional model for each species was generated using MAXENT 3.3.3 k  and included 75 % of species records for training and 25 % for testing the model. We employed receiver operating characteristic (ROC) analysis to measure model performance . The accuracy of each model prediction was tested by the area under the ROC curve (AUC) . AUC scores range from 0 to 1, with a score greater than 0.7 considered to be a good model performance .
Hair density on the abaxial leaf surface (HBL), hair density on the adaxial leaf surface (HDL), and hair density on the corolla outside (HCO) were highly correlated with each other (r > 0.9), and thus HBL and HCO were eliminated from subsequent multivariate analyses. Three vegetative characteristics, maximum leaf blade length (LL), maximum leaf blade width (LW), and hair density on the adaxial leaf surface (HDL), and eight floral characteristics, pedicel width (PW), corolla diameter (CLD), hair density on the corolla inside (HCI), corona diameter (CND), gynostegium length (GYL), gynostegium diameter (GYD), corpusculum length (COL), and pollinia length (POL) that exhibited significant differences (P < 0.05) among morphologically ‘pure’ V. atratum and V. japonicum populations (Table 1) were used for multivariate analyses. Inter- and intra-population variations in these morphological characteristics are shown in box plots (Additional file 1: Figure S2) and bar plots (Additional file 1: Figure S3). The intra-population variation range for these characteristics varied among putative hybrid populations. For example, PW values of the H01 population tended to be similar to those of V. atratum, whereas the PW values of the H02–H05 populations tended to be similar to those of V. japonicum. The ranges of variation in PW in the H06 and H07 populations were intermediate between those of V. atratum and V. japonicum.
The minimum, maximum, and median values of pollen fertility are summarized in Additional file 1: Table S5. All the populations examined exhibited high pollen fertility. Pollen fertility ranged from 97–100 % in V. atratum populations (A01–A03), 96–100 % in hybrid populations (H01–H07), and 100 % in V. japonicum populations (J01–J03).
A total of 144 alleles were found for 10 nrSSR loci in V. atratum populations compared with 165 alleles for V. japonicum. The mean number of alleles per locus (N a) within a population was 5.30–8.10 in V. atratum and 2.60–8.90 in V. japonicum (Additional file 1: Table S1). The allele richness and the expected heterozygosity were 3.89–4.82 and 0.597–0.697, respectively, in V. atratum, and 2.06–4.34 and 0.289–0.667, respectively, in V. japonicum (Additional file 1: Table S1). While most of the loci were in Hardy-Weinberg equilibrium (HWE), 31 of 260 cases showed significant deviations from HWE expectation (Additional file 1: Table S6). While most of the locus pairs did not show significant deviation from linkage equilibrium, one pair of loci in the H06 and J02 populations was in LD after applying the Bonferroni correction (P > 0.01) (Additional file 1: Table S7). All pair wise F ST estimates were significantly different from zero (P < 0.05). Genetic differentiation between all pairs of V. atratum and V. japonicum populations were moderate or large with F ST values ranging from 0.18 to 0.52 (Additional file 1: Table S8).
Population genetic structure
Long-term estimates of gene flow
Coalescent-based estimates obtained from the LAMARC program indicated that 4 Nm ‘pure’V. atratum →‘pure’ V. japonicum (19.27) was much higher than 4 Nm ‘pure’ V. japonicum→‘pure’V. atratum (1.18) (Additional file 1: Table S10). In addition, there was no evidence of unequal amount genetic input into the putative hybrid population, because of overlapping CIs for estimates of 4 Nm in both directions.
Species distribution modeling
Shifts in the ranges of V. atratum and V. japonicum with climate change
The potential ranges of V. atratum and V. japonicum predicted by SDM suggested that these two ecologically distinct species responded to climate changes in different ways across the Japanese Archipelago (Fig. 6). During the LGM, V. atratum was distributed mainly in the entire western part of Japan and along the Japan Sea side of Honshu, whereas V. japonicum was restricted to possible southern refugia, mainly at the southern tip of Kyushu and in a narrow coastal belt of the Pacific side of Honshu. During the warmer post-glacial period, V. atratum expanded its range across most of the territory of Honshu, whereas V. japonicum expanded its range across western Japan and some parts of eastern Japan.
Consequences of secondary contact and patterns of hybridization between V. atratum and V. japonicum
Morphological data confirmed the putative hybridization between V. atratum and V. japonicum. Morphologically ‘pure’ V. atratum and V. japonicum were readily discriminated based on the PCO analysis (Fig. 3) and variations in composite traits (Additional file 1: Figure S2–S3). All putative hybrid populations exhibited intermediate and/or parental-like phenotypes; however, morphological investigation alone is insufficient to explain gene flow-mediated hybridization patterns [63, 64]. We also used microsatellite markers to examine the genetic structure of hybrid populations. The Bayesian-based clustering analysis of nuclear microsatellite markers using STRUCTURE showed that the optimal value of K was 2 corresponding to the two morphologically ‘pure’ V. atratum populations and part of the morphologically ‘pure’ populations of V. japonicum (J01–J08) (Fig. 4). In the putative hybrid populations, most individuals exhibited V. atratum-like genotypes (H06 and H07), V. japonicum-like genotypes H02–H05 and admixed genotypes (H01). The DAPC of all individuals again showed two genetically distinct groups as did analysis based on STRUCTURE (Fig. 5). The H01, H06, and H07 populations were genetically closer to V. atratrum, whereas the H02–H05 populations were genetically closer to V. japonicum.
Most hybrid populations (H02–H05 and H06–H07) consist with parental-like genotypes. These genotypes may be generated by frequent backcrossing to one of parental species, because highly backcrossed individuals are difficult to distinguish from parental individuals due to the relatively small number of loci . Unidirectional backcrossing might be explained by reproductive system [8–11] and/or relative species abundance . Following reasons are possible to explain the fact. First, factors driving unidirectional backcrossing in plants usually include differences in flower structure , flowering time  and mating system  among interfertile species, differences in the offspring fitness of reciprocal crosses . In the present study, there is no difference in flower structure, flowering time and mating system between V. atratum and V. japonicum. In addition, biases in offspring fitness unlikely exist, because both V. atratum-like genotypes (H06 and H07), V. japonicum-like genotypes H02–H05 and admixed genotypes (H01) were detected. Alternatively, the unidirectional backcrossing might result from ‘pollen swamping’ by more abundant local species over later-colonizing species [66–68]. The colonizing species may have been subject to substantial pollen flow from the local species and repeated backcrosses with incoming pollen would generate abundant backcrossing to the local aboriginal species. Genetic constitutions of hybrid populations H06 and H07 can be explained by the prediction of colonization process (Fig. 6); during the LGM, V. atratum were already distributed in the area where the current H06 and H07 populations are located. After the glacial period, V. japonicum may have colonized this area from its southern refugia and been subject to substantial pollen flow from V. atratum. If backcrosses to V. atratum are selectively favored and do not exhibit reduced fitness relative to V. japonicum and V. atratum, they may ultimately displace ‘pure’ populations of V. japonicum and V. atratum.
Although the SDM results indicated that V. atratum colonized the area where the current H01–H05 populations were located earlier than V. japonicum (Fig. 6), the genetic constitutions of the H01–H05 populations differed from those of the H06 and H07 populations (Fig. 5). We consider that the genetic constitutions of these hybrid populations between two Vincetoxicum were influenced not only by the historical colonization process, but were also regulated by other factor, such as, abiotic environmental regimes. Variation in fine-scale habitats might promote differences in the strengths of selection on recombination genotypes and lead to variable spatial structure in a narrow region [15–22]. In the present study, population H01 show different genetic constitution with population H02-H05. During field investigation, we found hybrid populations H01 located in high cliff, whereas population H02-H05 located in rocky beaches. Thus, population H01 and H02–H05 may have suffered different ecological selection and show different genetic constitution.
Long-term hybrid swarm persistence
Hybrid swarms are believed to persist through two mechanisms [5, 69]. They may be maintained by frequent and strong gene flow that counterbalances selection against the hybrids. Alternatively, they may be self-sustaining and maintained largely by selective advantage. The hybridization between the two Vincetoxicum species undoubtedly falls into the latter circumstance. First, LD was observed in almost no combinations of loci in the hybrid populations (Additional file 1: Table S7) suggesting that the hybrid populations are not maintained by ongoing genetic input from parental populations and may have persisted for several generations. These attributes are characteristic of a hybrid swarm where the hybrids are self-sustaining . Second, the pollen fertility of all individuals in the hybrid populations was very high (Additional file 1: Table S5) providing high viability to facilitate mating between the individuals. Combined with the SDM results (Fig. 6), which predicted probable secondary contact during the post-glacial period, one possible scenario to account for the existence of the hybrid swarm is that the hybrid populations probably formed when V. atratum and V. japonicum came into contact during the post-glacial period and may be self-sustaining and have persisted for a long term – since V. atratum and V. japonicum came into contact during the post-glacial period.
The Bayesian-based clustering analyses showed signs of introgressive hybridization (Fig. 5). In the J09–J11populations, the morphologically ‘pure’ V. japonicum population was composed of genetic clusters of both V. atratum and V. japonicum. Coalescent-based estimates from LAMARC also indicated high historic gene flow among ‘pure’ V. atratum and ‘pure’ V. japonicum populations (Additional file 1: Table S10). The modeled distribution shows that during LGM, suitable habitat for V. japonicum J09–J11 populations still persist in where they exist currently. There would have been a situation that small V. japonicum populations continued to persist in the area (where J09-J11 are located currently) surrounded by much larger populations of V. atratum, resulting in directional introgression from the latter species.
Discrepancy in morphological and genetic identification of hybridity
The morphologies of individuals in the H06 and H07 hybrid populations were intermediate between those of V. atratum and V. japonicum (Fig. 3 and Additional file 1: Figure S2). Surprisingly, almost no individuals exhibited intermediate genotypes based on the results of microsatellite analyses (Fig. 4). A probable cause of this discrepancy is that morphology may be controlled by numerous genes, whereas a relatively small number of microsatellite loci were employed in the Bayesian-based assignment approach. Thus, microsatellite genotypes may be rapidly fixed to one of the parental species genotypes by backcrossing. However, the genes affecting morphology may have persisted over longer periods of time.
In this study, we demonstrated that climatic oscillations during the Quaternary period likely led to species range shift and subsequently secondary contact. We confirmed the hybridization between V. atratum and V. japonicum. Hybrid populations may be self-sustaining and have persisted since V. atratum and V. japonicum came into contact during the post-glacial period. Moreover, our findings show that the pattern of hybridization between V. atratum and V. japonicum is unidirectional and differs among populations. In addition, Bayesian-based clustering analysis and coalescent-based estimates of long-term gene flow showed patterns of introgressive hybridization in three morphologically ‘pure’ V. japonicum populations. We concluded that these differences in the genetic consequences of secondary contact are caused by historical colonization processes and/or natural selection.
Availability of supporting data
The datasets supporting the results of this article are available in the Dryad (doi:10.5061/dryad.mg5rs).
We thank Dr. S. Horie for technical advice and Mr. Y. Takahashi for sampling. We also thank the anonymous reviewers for helpful comments on the manuscript. This work was conducted under the framework of the “Precise Impact Assessments on Climate Change” of the Program for Risk Information on Climate Change (SOUSEI Program) supported by the Ministry of Education, Culture, Sports, Science, and Technology-Japan (MEXT). This study was also supported by a Grant-in-Aid to MM from the Japan Society for the Promotion of Science (JSPS KAKENHI 21570086, http://www.jsps.go.jp) and the Tohoku University Global Centre of Excellence (GCOE) Program “Center for Ecosystem Management Adapting to Global Change” (J03).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Comes HP, Kadereit JW. The effect of quaternary climatic changes on plant distribution and evolution. Trends Plant Sci. 1998;3(11):432–8.View ArticleGoogle Scholar
- Davis MB, Shaw RG. Range shifts and adaptive responses to Quaternary climate change. Science. 2001;292(5517):673–9.View ArticlePubMedGoogle Scholar
- Taberlet P, Fumagalli L, Wust-Saucy A, Cosson J. Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998;7(4):453–64.View ArticlePubMedGoogle Scholar
- Hewitt GM. Speciation, hybrid zones and phylogeography - or seeing genes in space and time. Mol Ecol. 2001;10(3):537–49.View ArticlePubMedGoogle Scholar
- Arnold ML. Natural hybridization and evolution. New York: Oxford University Press; 1997.Google Scholar
- Coyne JA, Orr HA. Speciation. 37th ed. Sunderland: Sinauer Associates; 2004.Google Scholar
- Petit RJ, Excoffier L. Gene flow and species delimitation. Trends Ecol Evol. 2009;24(7):386–93.View ArticlePubMedGoogle Scholar
- Field DL, Ayre DJ, Whelan RJ, Young AG. Patterns of hybridization and asymmetrical gene flow in hybrid zones of the rare Eucalyptus aggregata and common E. rubida. Heredity. 2011;106(5):841–53.PubMed CentralView ArticlePubMedGoogle Scholar
- Milne RI, Abbott RJ. Reproductive isolation among two interfertile Rhododendron species: low frequency of post-F1 hybrid genotypes in alpine hybrid zones. Mol Ecol. 2008;17(4):1108–21.View ArticlePubMedGoogle Scholar
- Wallace LE, Culley TM, Weller SG, Sakai AK, Kuenzi A, Roy T, et al. Asymmetrical gene flow in a hybrid zone of Hawaiian Schiedea (Caryophyllaceae) species with contrasting mating systems. PLoS One. 2011;6(9):e24845.PubMed CentralView ArticlePubMedGoogle Scholar
- Natalis LC, Wesselingh RA. Post-pollination barriers and their role in asymmetric hybridization in Rhinanthus (Orobanchaceae). Am J Bot. 2012;99(11):1847–56.View ArticlePubMedGoogle Scholar
- Lepais O, Petit RJ, Guichoux E, Lavabre JE, Alberto F, Kremer A, et al. Species relative abundance and direction of introgression in oaks. Mol Ecol. 2009;18(10):2228–42.View ArticlePubMedGoogle Scholar
- Currat M, Ruedi M, Petit RJ, Excoffier L. The hidden side of invasions: massive introgression by local genes. Evolution. 2008;62(8):1908–20.PubMedGoogle Scholar
- Petit RJ, Bodénès C, Ducousso A, Roussel G, Kremer A. Hybridization as a mechanism of invasion in oaks. New Phytol. 2004;161(1):151–64.View ArticleGoogle Scholar
- Cruzan MB, Arnold ML. Ecological and genetic associations in an Iris hybrid zone. Evolution. 1993;47(5):1432–45.View ArticleGoogle Scholar
- Hatfield T, Schluter D. Ecological Speciation in Sticklebacks: Environment-Dependent Hybrid Fitness. Evolution. 1999;53(3):866–73.View ArticleGoogle Scholar
- Hochwender CG, Fritz RS. Fluctuating asymmetry in a Salix hybrid system: The importance of genetic versus environmental causes. Evolution. 1999;53(2):408–16.View ArticleGoogle Scholar
- Johansen-Morris AD, Latta RG. Genotype by environment interactions for fitness in hybrid genotypes of Avena barbata. Evolution. 2008;62(3):573–85.View ArticlePubMedGoogle Scholar
- Johnston JA, Grise DJ, Donovan LA, Arnold ML. Environment-dependent performance and fitness of Iris brevicaulis, I. fulva (Iridaceae), and hybrids. Am J Bot. 2001;88(5):933–8.View ArticlePubMedGoogle Scholar
- Karrenberg S, Favre A. Genetic and ecological differentiation in the hybridizing campions Silene dioica and S. latifolia. Evolution. 2008;62(4):763–73.View ArticlePubMedGoogle Scholar
- Rieseberg LH, Sinervo B, Linder CR, Ungerer MC, Arias DM. Role of gene interactions in hybrid speciation: evidence from ancient and experimental hybrids. Science. 1996;272(5652):741–4.View ArticlePubMedGoogle Scholar
- Ross RIC, Agren JA, Pannell JR. Exogenous selection shapes germination behaviour and seedling traits of populations at different altitudes in a Senecio hybrid zone. Ann Bot-London. 2012;110(7SI):1439–47.View ArticleGoogle Scholar
- Ono Y, Igarashi Y. Natural history of Hokkaido. Hokkaido: Hokkaido University Press; 1991.Google Scholar
- Fujii N, Tomaru N, Okuyama K, Koike T, Mikami T, Ueda K. Chloroplast DNA phylogeography of Fagus crenata (Fagaceae) in Japan. Plant Syst Evol. 2002;232(1-2):21–33.View ArticleGoogle Scholar
- Sugahara K, Kaneko Y, Ito S, Yamanaka K, Sakio H, Hoshizaki K, et al. Phylogeography of Japanese horse chestnut (Aesculus turbinata) in the Japanese Archipelago based on chloroplast DNA haplotypes. J Plant Res. 2011;124(1):75–83.View ArticlePubMedGoogle Scholar
- Tono A, Iwasaki T, Seo A, Murakami N. Environmental factors contribute to the formation and maintenance of the contact zone observed in deciduous broad-leaved tree species in Japan. J Plant Res. 2015;128(4):535–51.View ArticlePubMedGoogle Scholar
- Qi X, Yuan N, Comes H, Sakaguchi S, Qiu Y. A strong ‘filter’ effect of the East China Sea land bridge for East Asia’s temperate plant species: inferences from molecular phylogeography and ecological niche modelling of Platycrater arguta (Hydrangeaceae). BMC Evol Biol. 2014;14(1):41.PubMed CentralView ArticlePubMedGoogle Scholar
- Mimura M, Mishima M, Lascoux M, Yahara T. Range shift and introgression of the rear and leading populations in two ecologically distinct Rubus species. BMC Evol Biol. 2014;14:209.PubMed CentralView ArticleGoogle Scholar
- Toyama H, Yahara T. Comparative phylogeography of two closely related Viola species occurring in contrasting habitats in the Japanese archipelago. J Plant Res. 2009;122(4):389–401.View ArticlePubMedGoogle Scholar
- Lee S, Maki M. Comparative phylogeographic study of Hosta sieboldiana and Hosta albomarginata (Asparagaceae) in Japan. Ecol Evol. 2013;3(14):4767–85.PubMed CentralView ArticlePubMedGoogle Scholar
- Saeki I, Dick CW, Barnes BV, Murakami N. Comparative phylogeography of red maple (Acer rubrum L.) and silver maple (Acer saccharinum L.): impacts of habitat specialization, hybridization and glacial history. J Biogeogr. 2011;38(5):992–1005.View ArticleGoogle Scholar
- Yamazaki T. Asclepiadaceae. In: Iwatsuki K, Yamazaki T, Boufford D, Ohba H, editors. Flora of Japan. 3ath ed. Tokyo: Kodansha; 1993. p. 168–82.Google Scholar
- Yamashiro T, Yamashiro A, Yokoyama J, Maki M. Morphological aspects and phylogenetic analyses of pollination systems in the Tylophora–Vincetoxicum complex (Apocynaceae-Asclepiadoideae) in Japan. Biol J Linn Soc. 2008;93(2):325–41.View ArticleGoogle Scholar
- Yamashiro T, Fukuda T, Yokoyama J, Maki M. Molecular phylogeny of Vincetoxicum (Apocynaceae-Asclepiadoideae) based on the nucleotide sequences of cpDNA and nrDNA. Mol Phylogenet Evol. 2004;31(2):689–700.View ArticlePubMedGoogle Scholar
- Chakraborty R, Weiss KM. Admixture as a tool for finding linked genes and detecting that difference from allelic association between loci. Proc Natl Acad Sci U S A. 1988;85(23):9119–23.PubMed CentralView ArticlePubMedGoogle Scholar
- Briscoe D, Stephens JC, O’Brien SJ. Linkage disequilibrium in admixed populations: applications in gene mapping. J Hered. 1994;85(1):59–63.PubMedGoogle Scholar
- Goodman SJ, Barton NH, Swanson G, Abernethy K, Pemberton JM, Goodman SJ, et al. Introgression through rare hybridization: A genetic study of a hybrid zone between red and sika deer (genus Cervus) in Argyll, Scotland. Genetics. 1999;152(1):355–71.PubMed CentralPubMedGoogle Scholar
- Maki M, Horie S, Yokoyama J. Comparison of genetic diversity between narrowly endemic shrub Menziesia goyozanensis and its widespread congener M. pentandra (Ericaceae). Conserv Genet. 2002;3(4):421–5.View ArticleGoogle Scholar
- Tada F, Yamashiro T, Maki M. Development of microsatellite markers for the endangered grassland perennial herb Vincetoxicum atratum (Apocynaceae–Asclepiadoideae). Conserv Genet. 2009;10(4):1057–9.View ArticleGoogle Scholar
- Nakahama N, Kaneko S, Hayano A, Isagi Y, Inoue-Murayama M, Tominaga T. Development of microsatellite markers for the endangered grassland species Vincetoxicum pycnostelma (Apocynaceae) by using next-generation sequencing technology. Conserv Genet Resour. 2012;4(3):669–71.View ArticleGoogle Scholar
- R Development Core Team. R: A Language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2012.Google Scholar
- Li Y, Itoi T, Takahashi H, Maki M. Morphological and genetic variation in populations in a hybrid zone between Leucosceptrum japonicum and L. stellipilum (Lamiaceae) in the central Japanese mainland. Plant Syst Evol. 2015;301(3):1029–41.View ArticleGoogle Scholar
- Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research--an update. Bioinformatics. 2012;28(19):2537–9.PubMed CentralView ArticlePubMedGoogle Scholar
- Kalinowski ST. Hp‐rare 1.0: A computer program for performing rarefaction on measures of allelic richness. Mol Ecol Notes. 2005;5(1):187–9.View ArticleGoogle Scholar
- Weir BS, Cockerham CC. Estimating F-Statistics for the Analysis of Population Structure. Evolution. 1984;38(6):1358–70.View ArticleGoogle Scholar
- Excoffier L, Lischer H. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10(3):564–7.View ArticlePubMedGoogle Scholar
- Rousset F. genepop’007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol Resour. 2008;8(1):103–6.View ArticlePubMedGoogle Scholar
- Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59.PubMed CentralPubMedGoogle Scholar
- Chen C, Durand E, Forbes F, François O. Bayesian clustering algorithms ascertaining spatial population structure: a new computer program and a comparison study. Mol Ecol Notes. 2007;7(5):747–56.View ArticleGoogle Scholar
- Earl DA. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4(2):359–61.View ArticleGoogle Scholar
- Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23(14):1801–6.View ArticlePubMedGoogle Scholar
- Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4(1):137–8.View ArticleGoogle Scholar
- Jombart T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24(11):1403–5.View ArticlePubMedGoogle Scholar
- Kuhner MK. LAMARC 2.0: maximum likelihood and Bayesian estimation of population parameters. Bioinformatics. 2006;22(6):768–70.View ArticlePubMedGoogle Scholar
- Kimura M, Ohta T. Stepwise mutation model and distribution of allelic frequencies in a finite population. Proc Natl Acad Sci U S A. 1978;75(6):2868–72.PubMed CentralView ArticlePubMedGoogle Scholar
- Rambaut A, Suchard MA, Xie D, Drummond AJ: Tracer v1.6. Available at http://beast.bio.ed.ac.uk/Tracer. 2014.
- Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25(15):1965–78.View ArticleGoogle Scholar
- Hijmans RJ, Guarino L, Cruz M, Rojas E. Computer tools for spatial analysis of plant genetic resources data: 1. DIVA-GIS. Plant Genet Resour Newsl. 2001;127:15–9.Google Scholar
- Phillips SJ, Dudík M. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography. 2008;31(2):161–75.View ArticleGoogle Scholar
- Fawcett T. An introduction to ROC analysis. Pattern Recogn Lett. 2006;27(8):861–74.View ArticleGoogle Scholar
- Fielding AH, Bell JF. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ Conserv. 1997;24(01):38–49.View ArticleGoogle Scholar
- Swets JA. Measuring the accuracy of diagnostic systems. Science. 1988;240(4857):1285–93.View ArticlePubMedGoogle Scholar
- Hardig TM, Brunsfeld SJ, Fritz RS, Morgan M, Orians CM. Morphological and molecular evidence for hybridization and introgression in a willow (Salix) hybrid zone. Mol Ecol. 2000;9(1):9–24.View ArticlePubMedGoogle Scholar
- Tsukaya H, Fukuda T, Yokoyama J. Hybridization and introgression between Callicarpa japonica and C. mollis (Verbenaceae) in central Japan, as inferred from nuclear and chloroplast DNA sequences. Mol Ecol. 2003;12(11):3003–11.View ArticlePubMedGoogle Scholar
- Boecklen WJ, Howard DJ. Genetic Analysis of Hybrid Zones: Numbers of Markers and Power of Resolution. Ecology. 1997;78(8):2611–6.View ArticleGoogle Scholar
- Ellstrand NC, Elam DR. Population Genetic Consequences of Small Population Size: Implications for Plant Conservation. Annu Rev Ecol Syst. 1993;24(1):217–42.View ArticleGoogle Scholar
- Levin DA, Francisco-Ortega J, Jansen RK. Hybridization and the Extinction of Rare Plant Species. Conserv Biol. 1996;10(1):10–6.View ArticleGoogle Scholar
- Rhymer JM, Simberloff D. Extinction by Hybridization and Introgression. Annu Rev Ecol Syst. 1996;27:83–109.View ArticleGoogle Scholar
- Anderson E. Hybridization of the Habitat. Evolution. 1948;2(1):1–9.View ArticleGoogle Scholar
- Latch EK, Kierepka EM, Heffelfinger JR, Rhodge OE. Hybrid swarm between divergent lineages of mule deer (Odocoileus hemionus). Mol Ecol. 2011;20(24):5265–79.View ArticlePubMedGoogle Scholar