Different Roles of Introgression on the Demographic Change in Two Snakebark Maples, Acer caudatifolium and A. morrisonense, with Contrasted Postglacial Expansion Routes

Hybridization frequently occurs in plant species. With repeated backcross, the introgression may influence evolutionary trajectories through the entry of foreign genes. However, the genetic admixture via hybridization events is often confused with the ancestral polymorphism, especially in closely related species that have experienced similar evolutionary events. In Taiwan, two independent-originated endemic snakebark maples have contrasted postglacial range expansion routes: northward and upward expansion in Acer caudatifolium and downward expansion in A. morrisonense. The range expansion causes the current parapatric distribution, increasing the possibility of introgression. This study elucidates how their genetic variation reflects introgression and historical demography. With 17 EST-SSR markers among the intensely sampled 657 individuals, we confirmed that the genetic admixture between species mainly was attributed to recent introgression instead of common ancestral polymorphism. The secondary contact scenario inferred by approximate Bayesian computation suggested that A. morrisonense received more genetic variations from A. caudatifolium. Introgression occurred in colonized Taiwan around the early Last Glacial Period. Furthermore, the demography of A. caudatifolium was more severely affected by introgression than A. morrisonense, especially in the wavefront populations with high altitude range expansion, implying an altitude-related adaptive introgression. In contrast, A. morrisonense exhibited ubiquitous introgression independent of postglacial expansion, suggesting that introgression in A. morrisonense was neutral. In terms of different genetic consequences, introgression had different demographic impacts on species with different altitude expansion directions even under the same climate-change conditions within an island.


Introduction
Hybridization in vascular plants is common [1] but not universal [2]. Perennials and outbreeding with characters of stabilized hybridity, such as agamospermy and vegetative spread, characterize hybrid hotspots [2]. Via hybridization, genes are invasive from other species [3], which is likely to prevent fixation [4] and slow the evolution of reciprocal monophyly [5] but less likely to cause the fusion of species [6]. Introgression can occur via repeated hybridization and backcross or through the easily permeable genome regions of species boundary. The former describes the phenomenon of incorporating alleles from one species into another, while the latter may, but not necessarily, imply the selectivity of the porous genome [5]. Heterogeneous recombination rates across the genome could also result in the semipermeable nature of species boundary [7]. In other words, introgression refers to the common phenomenon of foreign genes entering the gene pool of a species, but its mechanism may be complicated.     LGM) and the present, respectively. (c,d) denote the altitudinal distribution in the LGM and the present, respectively. The grid numbers of (c,d) were counted from (a,b), respectively. These plots showed the allopatric distribution of two maples in the LGM and the current parapatric distributions. Data on ENM distributions are adopted from Luo et al. [22] and Chang et al. [23].

Sampling and DNA Extraction
In total, we sampled 657 samples in Taiwan for SSR genotyping, in which 371 individuals of A. caudatifolium were collected from 20 populations, and 286 individuals of A. morrisonense were from 19 populations (Table S1). The genomic DNA was extracted from the fresh leaves using modified Cetyltrimethylammonium bromide (CTAB) approach [24,25]. Polyvinylpyrrolidone (PVP) and Polyvinylpolypyrrolidone (PVPP) were added to remove polysaccharides, phenolic compounds, and other secondary metabolites. RNaseA was used to remove excess RNA. NanoDrop microvolume spectrophotometer was used to qualify and quantify the extracted DNA.

Transferable Expressed Sequence Tag-Simple Sequence Repeat (EST-SSR)
The primers for the microsatellite DNA (simple sequence repeats, SSR) were designed from putative transcripts (Source Version: 01052015) of A. saccharum in Hardwood Genomics Project (https://www.hardwoodgenomics.org/) (accessed on 24 July 2017) using SciRoKo v3.4 [26]. We set the criteria of two-or-more nucleotide motifs with at least four repeats to select the SSR loci. All selected loci were compared with the transcriptome of A. negundo (Accession: VFFP) in One Thousand Plants (1KP) Consortium (https://sites.google.com/a/ualberta.ca/onekp/) (accessed on 30 August 2017), and the E-value must be <1 × 10 −10 . The lengths of all SSR loci were <500 bps with 40~60% GC content. The designed primers were checked with IDT Oligo Analyzer 3.1 (http://sg.idtdna.com/calc/analyzer) (accessed on 20 September 2017), and the differences of the annealing temperatures were <3 • C between forward and reverse primers. In total, we designed 47 primer pairs, of which 29 primer pairs were amplifiable. However, nine of 29 SSR loci were monomorphic between species, and three of the remaining generated >70% null alleles. A total of 17 available SSR loci were obtained finally. All primer pairs were published in Chang et al. [23].
PCR products of these 17 loci were genotyped by capillary electrophoresis on an Applied Biosystems 3730 DNA Analyzer (Applied Biosystem, San Francisco, CA, USA). Fragment sizes of SSR alleles were read by Peak Scanner version 1.0 (Applied Biosystem, San Francisco, CA, USA) and corrected by reference to the size standard ABI GS500 LIZ (Applied Biosystems, San Francisco, CA, USA) with the assistance of the National Center for Genome Medicine, Academia Sinica, Taiwan. Peaks with a minimum height of 100 were read, and the sizes falling into the expected range were manually checked and adjusted. Genotypes of A. morrisonense were used for the demographic analyses and published in Chang et al. [23].

Genetic Diversity Estimation
The observed (Ho) and expected heterozygosity (He), Wright's inbreeding coefficient (F IS ) of the SSR loci of each species were calculated by GenAlEx 6.502 [27]. The private alleles of each species were calculated. Linear regression was conducted to test the geographic relationship (altitude, longitude, and latitude) with private allele frequency.

Testing Introgression by STRUCTURE
The genetic structure between the two species was assessed using the Bayesian clustering analysis (BCA) by STRUCTURE v. 2.3.4 [28,29]. The grouping numbers K = 1-21 were set with ten replications. Each replication conducted a 10 6 -times Markov chain Monte Carlo simulation with 20% burn-in. The best K was evaluated by ∆K and CV error in Structure Harvester [30]. We defined the admixture coefficient (Q) > 0.05 as the introgression.

Approximate Bayesian Computation (ABC)
Four evolutionary scenarios addressing the species divergence and interspecific introgression between parapatric species referred to [31] were tested by Approximate Bayesian Computation (ABC) using ABCtoolbox [32]: (1) strict isolation (SI) model: no gene flow  (Figure 3d). These four scenarios described all possible situations in which two lineages coalesced back into a common ancestor and are mutually exclusive. The parameter setting of models is shown in Table S2. We performed 10 6 simulations in fastsimcoal2 [33], and the best 2000 simulations were retained using the R package abc [34] to perform model selection with the neuralnet method. The model with the highest posterior probability (PP) and Bayes factors was suggested as the optimal scenario. Goodness-of-fit test was also used to confirm the preferred model for the observed data. Parameters of effective population size (Ne anc , Ne cau , Ne mor : effective population sizes of the common ancestor, A. caudatifolium, and A. morrisonense, respectively), migration rate (m 1 , m 2 ), and divergence time (t 1 ), and the time to gene flow (t 2 ) of the optimal model were calculated. On the basis of these four scenarios, we allowed a changing effective population size (N e ) after species divergence at t 2 . gression.

Approximate Bayesian Computation (ABC)
Four evolutionary scenarios addressing the species divergence and interspecific introgression between parapatric species referred to [31] were tested by Approximate Bayesian Computation (ABC) using ABCtoolbox [32]: (1) (Figure 3d). These four scenarios described all possible situations in which two lineages coalesced back into a common ancestor and are mutually exclusive. The parameter setting of models is shown in Table S2. We performed 10 6 simulations in fastsimcoal2 [33], and the best 2000 simulations were retained using the R package abc [34] to perform model selection with the neuralnet method. The model with the highest posterior probability (PP) and Bayes factors was suggested as the optimal scenario. Goodness-of-fit test was also used to confirm the preferred model for the observed data. Parameters of effective population size (Neanc, Necau, Nemor: effective population sizes of the common ancestor, A. caudatifolium, and A. morrisonense, respectively), migration rate (m1, m2), and divergence time (t1), and the time to gene flow (t2) of the optimal model were calculated. On the basis of these four scenarios, we allowed a changing effective population size (Ne) after species divergence at t2. The branch thickness represents the effective population size (Ne) in the recent (t2), early (t1), and of the most recent common ancestor (i.e., the uppermost part where the branches converge). Horizontal arrows indicate the direction of gene flows from Acer caudatifolium to A. morrisonense (m1) and the opposite direction (m2). Model parameters Ne and m denote the effective population size and migration rate. In these scenarios, two species diverged at t1 ago, and their population sizes were allowed to change at t2, the time of stopping interspecific gene flow in (AM) scenario or starting to gene flow in the (SC) scenario. Abbreviations of four scenarios: (SI): strict isolation, no interspecific gene flow was allowed after species divergence; (AM): ancient migration, interspecific gene flow was allowed at the early stage of species divergence; (SC): secondary contact, interspecific gene flow was allowed at recent periods; (CM): continuous migration, interspecific gene flow was allowed at all times after species divergence. The branch thickness represents the effective population size (N e ) in the recent (t 2 ), early (t 1 ), and of the most recent common ancestor (i.e., the uppermost part where the branches converge). Horizontal arrows indicate the direction of gene flows from Acer caudatifolium to A. morrisonense (m 1 ) and the opposite direction (m 2 ). Model parameters N e and m denote the effective population size and migration rate. In these scenarios, two species diverged at t 1 ago, and their population sizes were allowed to change at t 2 , the time of stopping interspecific gene flow in (AM) scenario or starting to gene flow in the (SC) scenario. Abbreviations of four scenarios: (SI): strict isolation, no interspecific gene flow was allowed after species divergence; (AM): ancient migration, interspecific gene flow was allowed at the early stage of species divergence; (SC): secondary contact, interspecific gene flow was allowed at recent periods; (CM): continuous migration, interspecific gene flow was allowed at all times after species divergence.

Genetic Diversity
Each of A. caudatifolium and A. morrisonense have one intraspecific monomorphic locus among the 17 interspecific polymorphic SSR loci. There is no significant difference in the genetic diversity between these two species estimated by SSR: the Ho and He are 0.169 ± 0.044 and 0.341 ± 0.073 in A. caudatifolium and 0.180 ± 0.017 and 0.261 ± 0.062 in A. morrisonense, respectively (Table 1). However, relatively abundant private alleles were detected in A. caudatifolium than A. morrisonense (Table 1)  sonense, respectively (Table 1). However, relatively abundant private alleles were detected in A. caudatifolium than A. morrisonense (Table 1). A further linear regression showed a positive association of private allele frequency of populations with altitude in A. caudatifolium (t = 2.699, p = 0.024, Figure 4a), but not in A. morrisonense (t = −0.372, p = 0.725, Figure  4a). Neither longitude nor latitude was related to the private alleles (Figure 4b,c). Detailed estimates of genetic diversity in each population are listed in Table S3. The pairwise FST among populations between the two species ranges 0.249-0.718. The average pairwise FST among populations of A. caudatifolium (FST = 0.144) is higher than A. morrisonense (FST = 0.106, p = 1.5 × 10 −11 , Table S4).  Analysis of molecular variance showed a significant divergence between two species (ΦCT = 0.478, p < 10 −5 ) and among populations (ΦSC = 0.149, p < 10 −5 ) ( Table 2). In terms of individual species, most genetic variation existed within populations (82.5% and 89.9% in A. caudatifolium and A. morrisonense, respectively), but there was still significant genetic differentiation between populations (ΦST = 0.175 and 0.101 in A. caudatifolium and A. morrisonense, respectively, Table 2). The significant population differentiation reflected the positive inbreeding coefficient estimates (F) in both species. Analysis of molecular variance showed a significant divergence between two species (Φ CT = 0.478, p < 10 −5 ) and among populations (Φ SC = 0.149, p < 10 −5 ) ( Table 2). In terms of individual species, most genetic variation existed within populations (82.5% and 89.9% in A. caudatifolium and A. morrisonense, respectively), but there was still significant genetic differentiation between populations (Φ ST = 0.175 and 0.101 in A. caudatifolium and A. morrisonense, respectively, Table 2). The significant population differentiation reflected the positive inbreeding coefficient estimates (F) in both species.

Genetic Introgression
The BCA showed an apparent genetic grouping between the two maples. However, some individuals still had a few genetic components of another species. Five individuals of A. morrisonense were suggested to be misidentifying in morphology because of the high genetic component of A. caudatifolium (92.29~99.80%), and one sample with 65.74% component of A. caudatifolium was assigned as a hybrid individual. Taking Q > 0.05 as the threshold, 9 of the 371 A. caudatifolium (2.43%) and 15 of the 281 A. morrisonense (5.34%) were genetically introgressive. Notably, these introgression events were not confined to specific populations but were occurred in sympatric (altitudinal parapatric) populations. This analysis implied that introgression was recent, common, and seemed asymmetric from A. caudatifolium toward A. morrisonense in Taiwan.
When mapping the genetic components along longitude and latitude, we found that introgression in both species spanned a wide range of longitudes (120.

Approximate Bayesian Computation
The genetic admixture that occurred early or recently after species divergence was tested in four evolutionary scenarios of isolation-with-migration by the approximate Bayesian computation (ABC) (Figure 3). Among the four scenarios, the strict isolation (SI) that did not allow any gene flow between species achieved the lowest posterior probabil-

Approximate Bayesian Computation
The genetic admixture that occurred early or recently after species divergence was tested in four evolutionary scenarios of isolation-with-migration by the approximate Bayesian computation (ABC) (Figure 3). Among the four scenarios, the strict isolation (SI) that did not allow any gene flow between species achieved the lowest posterior probability (PP = 0.005), while the secondary contact (SC) allowing interspecies gene flow at the late stage of divergence achieved the highest PP (PP = 0.488). The continuous migration (CM, PP = 0.360) and ancient migration (AM, PP = 0.148) scenarios achieved the second and third high PP, respectively. The extremely smaller PP in the SI scenario from the other three scenarios with gene flow showed that the gene flow continued after speciation, i.e., porous genomes between these two maples. The SC scenario also provides a good fit to the observed data (p = 0.45, Figure S1). Pairwise model comparisons between scenarios by Bayes factor were listed in Table S5.

Secondary Contact during Colonization to Taiwan in Last Glacial Period
Past genetic hybridization studies on Acer section Macrantha mainly focused on phylogenetics and hybrid ancestry [18]. This study, via the best ABC scenario "Secondary Contact" and introgression components inferred by BCA, evidenced that genetic admixture comes from recent and ongoing introgression between A. caudatifolium and A. morrisonense rather than the ancestral polymorphisms. This continued introgression is not surprising, as A. caudatifolium and A. morrisonense bloom synchronously from March to April within a confined space of an island. Since these two island maples were allopatrically distributed during the LGM until the middle-Holocene range overlapping [22,23], we speculate that the recent introgression is due to the secondary contact since postglacial expansion.
However, the onset time of introgression (t 2 ) was longer than the end of LGM. In South Asia, the sea level rose at the rate of 0.41 m/100 yrs between 19.0 kya to 14.6 kya and accelerated to 5.33 m/100 yrs between 14.6 and 14.3 kya [35], suggesting the end of LGM at approximately 14~19 kya in Asia. Under the SC scenario, t 2 was 18,870 generations ago, which must be earlier than the LGM unless there was one year per generation. However, although snakebark maples grow faster, according toĎurkovič [36], A. caudatifolium takes at least five years to mature and flower. In other words, A. caudatifolium and A. morrisonense started introgression much earlier (c. 94 kya at the beginning of the Late Pleistocene) than the end of LGM if considering five years per generation.
The onset time of introgression roughly matches the beginning of the Last Glacial Period (LGP) of the Quaternary (c. since 110 kya), during which the sea level fell > 50 m lower than the present day [37], resulting in the exposure of the seabed (average depth < 50 m at present) of the Taiwan Strait (i.e., the Fujian-Taiwan land bridge or, namely, Dōngshān land bridge). Although unsure of the exact time, the land bridge formation benefited the island colonization of two maple species. We speculate that the interspecific gene flow is related to entering Taiwan. In other words, this secondary contact inferred by ABC does not refer to the post-LGM range expansion, but to the events of colonization to Taiwan.

Rapid Range Expansion Rather Than Hybridization Shaped the Genetic Diversity of Wavefront Populations
The founder effect of colonizing Taiwan also explains the severe population decline from t 1 to t 2 inferred by ABC. Although previous studies only backtracked paleodistributions in LGM [22,23], these two species might also be confined within refugia since the earlier stage of the Quaternary LGP because of the continuous cold and aridity. However, we detected high frequencies of private alleles (Table 1), implying a rapid recovery of population size. We speculated that the reduced genetic diversity through population decline was supplemented by introgression, where the abundant private alleles reflect the rare-allele phenomenon after introgression [38]. If this hypothesis is supported, the population decline will be mitigated, viz., smaller Ne t1 /Ne t2 in the SC scenario (allowing introgression) than in the SI scenario (not allowing interbreeding).
However, the private allele frequencies were not high among hybrid populations (Table  S3). In addition, the mitigation of population decline was only detected in A. caudatifolium but not in A. morrisonense. Unlike the prevalent introgression in A. morrisonense, A. caudatifolium, with lower introgression rates, received foreign alleles only at higher altitudes. If introgression can complement the reduced genetic diversity, the relaxation of demographic constraint should be more pronounced in A. morrisonense. Therefore, the hypothesis of genetic rescue by introgression could be rejected, and the abundant private alleles were alternatively explained by the rapid range expansion after LGM. Especially for A. caudatifolium, the positive association of private allele with altitude reflects its post-LGM upward expansion [22] because of the higher mutation rate in the wavefront of expanding populations [39].
Nevertheless, the population decline of A. caudatifolium under the SC model was indeed lessened than that without introgression. We also found more genetic admixture at high elevation, the wavefront of population expansion. Spatial expansion surfs not only rare alleles into newly colonized populations but also accelerates native genes into the invading species [17]. In other words, introgression may advance the upward expansion by supplying new genetic variation in A. caudatifolium. High proportional new alleles at high elevation may also reflect the altitude-related adaptive introgression in A. caudatifolium (Figure 4a). The introgression and increasing rare alleles are conducive to widening niches of frontier populations. Such a phenomenon was not rare. For example, alpine Carex curvula adapts to edge environments in marginal habitats by introgressive hybridization [14]; the Mexican teosinte introgressed from maize expanded its niches to facilitate the invasion in Europe [40]. Introgression improves the fitness of postglacial expanding species within the space-constraint island.
However, an abundance of rare alleles and introgressions in wavefront populations was not observed in A. morrisonense. In the ABC simulations, its population size declined more severely in the SC scenario than in the SI scenario, suggesting that introgression did not rescue the irreparable genetic loss of critical range shrinkage in the middle Holocene [23]. High proportional private alleles were not detected in wavefront and introgression populations and were related to neither longitude, latitude, nor altitude. Thus, despite being more receptive to foreign genes than A. caudatifolium, it is not reflected in adaptiveness as with A. caudatifolium. Chang et al. [23] also suggested that adaptation-independent pollen flow homogenized genetic differences between distant populations whose seed colonization was restricted by mountain barriers. This also indicates that A. morrisonense is more receptive to pollens of foreign genotypes, implying a higher chance of neutral introgression.

Conclusions
In conclusion, this study proves that the genetic admixture of A. caudatifolium and A. morrisonense was attributed to the recent introgression of colonizing Taiwan since the LGP instead of from common ancestral polymorphism. However, the introgression did not recover the genetic loss of the LGP bottleneck effect. The increasing private alleles can only explain the post-LGM range expansion but not the rare allele phenomenon of introgression.
However, the high-elevational introgression may adaptively facilitate upward range expansion of A. caudatifolium, but such a phenomenon was not found in A. morrisonense whose introgression seems to be neutral and unrelated to its postglacial range expansion. This study outlined the spatiotemporal genetic changes of the introgression between these two island endemic maple trees following entering Taiwan and presented arguments for their different adaptive propensities during introgression.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11050644/s1, Figure S1: Histogram of the distribution under secondary contact (SC) scenario and the observed value; Figure S2: The frequency of marginal density of the best-fit scenario "secondary contact" (SC); Table S1: Sampling sites and sample size; Table S2: Search ranges of parameters in ABC model; Table S3: Genetic diversity of each sampled population and species;