Scutellaria barbata is distributed in northern Taiwan. Some populations were recorded from western and northeastern Taiwan. This species is usually scattered in the sallow hills, on farmland ridges and along roadsides, and is usually distributed sporadically, making it hard to form a dense population. Scutellaria taipeiensis is very rare and restricted to the Taipei Basin of northern Taiwan. It grows attached to bare rock or the surrounding soil, usually in sunny environments along trails (Huang et al., 2003). Therefore, each population size of these two species is quite limited and uneven through their distribution range. We collected samples from populations of both bar and tpe distributed in the Taipei Basin, and the surrounding areas. One population of bar was sampled in southwestern Taiwan to test the effect of long-distance geographic and climatic differences. In total, six and three populations of bar and tpe were collected, respectively (Table 1). Because of their morphological similarity, we first sampled in the MK population, which is the type location of tpe; the collected samples were subsequently checked and confirmed by one of the original authors of S. taipeiensis (Mr. Arthur Hsiao). Identification of other tpe populations was based on the macromorphology and growth habit of the MK population. In total, 80 and 74 individuals of bar and tpe were sampled, respectively, for genotyping and sequencing. The seed coat patterns of all populations were further checked by Scanning Electron Microscopy (SEM).

Seeds were collected from all studied populations for SEM observation. An accelerating voltage of 15 kV was applied in order to obtain high image resolution signals. SEM was conducted using the tabletop SEM TM3000 (Hitachi, Tokyo, Japan).

Chloroplast DNA fragments ndhF-rpl32 and rpl32-trnL were chosen for sequencing. The PCR amplification and sequencing conditions followed Chiang et al.'s (2012a) protocol. According to Chiang et al.'s (2012b) primer test, 14 microsatellite primers were amplifiable for both bar and tpe. Herein these 14 primer sets were used to test for polymorphism at the population level; three monomorphic loci were discarded in the further analyses. The amplification and genotyping conditions followed Chiang et al.'s (2012b) protocol. All sequences were deposited in GenBank (Supplementary Table 1).

After sequence alignment and trimming the 5′ and 3′ ends to equal lengths, neutrality tests were conducted by Tajima's (1989) D-test with coalescent simulations. For the neutrality test of microsatellite loci, we used the Fdist (Beaumont and Nichols, 1996) and Bayesian-based outlier analysis of the multinomial-Dirichlet model (BayeScan analysis) to assess the F_ST distribution among all loci. The genetic diversity of chloroplast DNA sequences was then estimated by DnaSP 5.10 (Librado and Rozas, 2009). The genetic diversity of microsatellite loci was estimated using Arlequin 3.5 (Excoffier and Lischer, 2010) and GenAlEx 6.5 (Peakall and Smouse, 2012).

The analysis of molecular variance (AMOVA) was used for testing the degrees of genetic variation between species, among, and within populations. We conducted the AMOVA using Arlequin 3.5 (Excoffier and Lischer, 2010).

Bayesian clustering analysis (BCA) was used to examine the similarity and divergence of genetic components among populations. The BCA was performed using STRUCTURE 2.3.4 (Hubisz et al., 2009). A non-admixture model was applied to STRUCTURE with a priori sample localities. The posterior probability of grouping number (K = 1–10) was estimated by 10 independent runs using one-million steps Markov chain Monte Carlo (MCMC) replicates after a 100,000-step burn-in for each run to evaluate consistency. The best grouping number was evaluated by ΔK (Evanno et al., 2005) in STRUCTURE HARVESTER v. 0.6.94 (Earl and Vonholdt, 2012).

The discriminant analysis of principal components (DAPC) was used for distinguishing bar and tpe and the populations that were subjectively identified. Properties of the “without a priori” using partial synthetic variables to minimize variation within groups (Manel et al., 2005; Jombart et al., 2010) could help to evaluate the artificial classification objectively. The DAPC was conducted using the R (R Core Team, 2015) package adegenet (Jombart, 2008). Kruskal–Wallis tests for the first two principal components (PCs) and the first two linear discriminants (LDs) of the DAPC were conducted to test the genetic divergence between species and between populations.

Based on the estimated genetic variation and genetic structure and current geographic distributions, two evolutionary scenarios were proposed: the first is the demographic bottleneck of both bar and tpe under common environmental pressures with recent interspecific gene flow (secondary contact) after completed speciation (Figure 2A); the second scenario illustrates the recent founder speciation of tpe derived from the progenitor species bar with continuous gene flow (i.e., incomplete speciation, Figure 2B). We further tested two extended scenarios of the recent founder speciation scenario, namely that tpe originated multiple times (Figure 2C) or in a single origin from bar (Figure 2D). Testing these two extended scenarios helps clarifying whether the species-determining morphological traits are derived from the common ancestor or are the consequence of convergence. An approximate Bayesian computation (ABC) was used to test these evolutionary scenarios. First, package simcoal2 of the ABCtoolbox was used to generate one million pseudodata for simulation (Wegmann et al., 2010). We used arlsumstat to estimate diversity parameters (e.g., number of alleles, allele frequency, fixation indices F_ST, F_IS, F_IT, etc.) of the pseudodata and transform them by the pls. To validate the best scenario, the general linear model (GLM) post-sampling regression adjustments on the best 5,000 simulated parameters were retained for obtaining the marginal density of models by the ABCestimator. The best evolutionary scenario was selected according to the Bayes factor and the ratio of marginal densities.

Grinnellian niches are defined as the scenopoetic environmental variables a species requires to survive, such as the temperature and precipitation, etc. (Grinnell, 1917). Nineteen bioclimatic variables (Bioclim) that represent the Grinnellian niches were extracted from the WorldClim website (http://www.worldclim.org/bioclim). Principal component analysis (PCA) of independent climatic variables to reduce the dimensionality defining the niche space allowed comparison of the integrity of Grinnellian niches between bar and tpe. Firstly, the variables with a high variance inflation factor (VIF) were removed to reduce the multicollinearity. This procedure was repeated until the VIF values of all remaining variables were < 10 and finally six bioclimatic variables were retained for the further analyses. Species occurrence data were collected from the Global Biodiversity Information Facility (GBIF, http://www.gbif.org), herbarium records, public records available on the web, and our collecting sites. Unclear records and wrong identifications were ruled out. Finally, 56 records of bar and five records of tpe were used for the PCA. We further used multinomial logistic regression (MLGR) to test the effect of bioclimatic variables on the prediction of species occurrence. Because the simpler model of no correlation among predictors (bioclimatic variables) could not be rejected by the model of correlated predictors (likelihood ratio statistic = 12.214, df = 37, P = 0.99996), we used the model with independent predictors for MLGR analysis. The significance of each predictor was tested by the type-II analysis of variance (ANOVA) Wald χ² test.

To predict the potential habitats of both bar and tpe, species distribution models were built under the maximum entropy model implemented in MaxEnt (Phillips and Dudík, 2008) and the R package dismo (Hijmans et al., 2013). The occurrence data of both species were obtained from our sampling sites, herbarium data, and GBIF. The bioclimatic variables were also downloaded from WorldClim with a resolution of 30 arc-sec (about 1 × 1 km). A maximum of 2,000 iterations of each species were conducted, and species occurrence data were randomly divided (10%) to train the model. One regularization multiplier and 10,000 background points were set to create models for each species set. The logistic output, consisting of a grid map with a suitability value range from 0 to 1, was generated and visualized with the R package maptools (Lewin-Koh et al., 2011).

To test how the geographic distance and environmental differences affect the genetic composition, the Mantel test of Spearman correlation was performed between genetic, geographic, and environmental (climatic) distances. Pairwise F_ST was calculated between populations based on Nei's (1977) approach and then Rousset's (1997) estimate F_ST/(1 − F_ST) was used as the genetic distance metric. Geographic distance was estimated using the Euclidean distance according to three dimensional factors (latitudes, longitudes, and altitudes). Environmental distance was also calculated using the Euclidean distance with six retained Grinnellian niches variables. A partial Mantel test between genetic and environmental distances controlled for the geographic distance was also done. A multiple matrix regression with randomization (MMRR) was further performed to test whether the genetic distance responds to the variation of geographic and/or environmental distances. We first examined the correlation between the genetic distance and geographic distance and between genetic distance and environmental distance. Autocorrelation between the geographic and environmental distances was further examined. Finally, the joint effect of both geographic and environmental distances on genetic distance was examined. Regression coefficients of Mantel test (ρ) and MMRR (β) and their significance were determined based on 9,999 random permutations.

To understand whether the climatic variables (Grinnellian niches) explain the genetic distribution of populations of bar and tpe, partial distance-based redundancy analyses (partial dbRDA) were performed. The first five principal components (PC1–PC5) were used to represent the genetic components, and six retained climatic variables were used as predictors conditioning on the geographic distribution (latitude and longitude). Analysis of variance was further used to test the significance of each predictor. The partial dbRDA was conducted using the R package vegan. We further analyzed the distribution of climatic variables along ordination axes using the GLM. We also used the MLGR to test the effect of bioclimatic variables on populations. The type-II ANOVA Wald χ² test was used to test the significance of each bioclimatic factor.

Because the characters of vegetative morphology that are the most distinctive between bar and tpe are variable and can be affected by growth and nutrient conditions, we specifically checked the seed coat patterns, another distinguishing trait in the reproductive organs, to confirm the stability of this character for species delimitation. Based on our observations, there is slight variation in the seed size among individuals but no obvious difference between species (Figure 3). The SEM also shows that the nutlet coats of mature seeds in all populations of bar and tpe are of the rounded concentric type (i.e., type 3 described in Huang et al., 2003, Figure 3). The radiated umbrella-like nutlet coat (type 2) described by Hsieh and Huang (1995) was only seen in immature seeds (the bottom-left picture of Figure 3). After repeatedly comparing the nutlet-coat patterns, we believe that the differences between these two character states (i.e., rounded-concentric vs. radiated umbrella-like) reflect differences in seed maturity rather than differences between species.

Before estimating genetic diversity, we performed neutrality tests by Tajima's D-test for cpDNA sequences and by Fdist and the multinomial-Dirichlet model (BayeScan) for the outlier analysis for the microsatellite loci. No deviation from zero in Tajima's D (D = 1.315, P > 0.10, coalescent simulation D = −0.066, 95% confidence interval −1.550–2.112, P = 0.911) and no positive or negative outliers of microsatellite loci (Supplementary Figure 2) suggest that the genetic markers used in this study all evolve neutrally and are appropriate for further genetic population assessment and speciation model tests.

In cpDNA, only the population YH of bar and populations MK and EG of tpe revealed nucleotide polymorphisms. The genetic diversity index Π of these three populations (YH, MK, and EG) is 0.143, 0.154, and 0.051, respectively, and the index θ_W is 0.314, 0.322, and 0.237, respectively. In total, when considering the indels six haplotypes were obtained, and the haplotype frequencies, distributions, and relationships (i.e., the minimum spanning tree) are illustrated in Figure 4. Genetic diversity assessed by microsatellite markers revealed varying degrees of genetic diversity among populations, ranging from 0.105 ± 0.039 to 0.296 ± 0.088 in expected heterozygosity (H_e) in bar and from 0.189 ± 0.080 to 0.230 ± 0.069 in H_e in tpe (Table 2). Among the 11 loci examined, roughly half are polymorphic (5.556 ± 1.878 polymorphic loci) in every population and the high proportions of private alleles per locus (ranges from 0 to 0.364 ± 0.203, Table 2) indicate a genetic fixation phenomenon in small populations by the drift effect. The drift effect on small populations is reflected in the significantly high genetic differentiation among populations (F statistic = 0.681, P < 0.00001) with a relatively small proportion of genetic variation within populations (35.94% variation, F = 0.641, P < 0.00001; Table 3). However, genetic variation does not contribute to species divergence, i.e., genetic differentiation between species is lacking (F < 0, P = 0.820, Table 3).

To understand the pattern of genetic distribution of these two species, we performed a BCA in STRUCTURE and ordination analysis by DAPC. The STRUCTURE result suggested the best grouping number (K) is 3 and the second-best K is 2 based on the ΔK (Figure 5A). When K = 2, population EG of tpe and population WL of bar were inferred to be composed of identical genetic components, and population YL of bar is composed of both genetic components (Figure 5B). When K = 3, DA of tpe and YL of bar were inferred to be two further components different from EG and WL, and the other populations were comprised of both genetic components of DA and YL (Figure 5C). These two inferences are similar and consistently indicate that the populations WL and EG are divergent from the other populations, and such divergence is inconsistent with the species delimitation.

In the DAPC analysis, which uses one discriminant function to distinguish five principal components (PCs), a broader range of genetic variation was detected in bar than in tpe, and the range of genetic variation of tpe is skewed but within the range of bar at the species level (Figure 6A). A further DAPC analysis on populations shows a similar result, i.e., that the range of genetic variation of three populations of tpe was clustered within the range of bar. Two main clusters were distinguished by two discriminant functions for five PCs: one is the cluster of the WL of bar and EG of tpe, the other populations belonging to another cluster (Figure 6B). The clustering pattern inferred by STRUCTURE and DAPC is congruent.

We further tested the differences in genetic variation according to the transformed variables of the first two principal components PC1 and PC2 of PCA and the first two linear discriminants LD1 and LD2 of DAPC at both the species level and the population level. Only the transformed variable PC2 (Kruskal–Wallis test, χ² = 32.33, df = 1, P = 1.301e–08) revealed species divergence between bar and tpe but no species-level differences were found in PC1 (P = 0.872), LD1 (P = 0.810), or LD2 (P = 0.864; Figure 7). At the population level, all four transformed vectors showed significant differences in genetic variation among populations (Kruskal–Wallis test, P < 2.2e–16, Figure 7). Taken together, these results suggest obvious population divergence but no or less genetic differentiation between species.

The non-significant genetic differentiation between species may suggest frequent gene flow geographically after secondary contact, or alternatively, the sharing of large amounts of ancestral polymorphism. Following these two possibilities, evolutionary scenarios of early divergence with secondary contact (Figure 2A) and recent founder speciation with incomplete divergence (Figure 2B) were tested using ABC. The former scenario (secondary contact) scored a marginal density of 9.24 × 10⁻¹³ (0.199%), much lower than the founder-speciation scenario (marginal density: 4.63 × 10⁻¹⁰, 99.801%). The estimated ancestral population size of bar is suggested to have decreased under the founder-speciation scenario (N_A = 597.487, 95% CI: 164.557–8172.980; N_B' = 199.497, 95% CI: 126.479–4701.65, Figure 8A). However, a recovery of the current population size was estimated (N_B = 683.417, 95% CI: 235.953–958.165), but the founder populations of tpe have a very small effective population size (N_T < 100, Figure 8B). The estimated time since the founder speciation of tpe is less than 100 generations (Figure 8D), indicating the poor inference of species divergence despite no or little interspecific gene flow (m < 0.01, Figure 8C) between bar and tpe. The estimated mutation rate (μ) of each locus per generation in secondary contact scenario has been estimated (μb = 6.965 × 10⁻⁴, 95% CI: 2.238-9.703 × 10⁻⁴, Figure 8H).

Since the divergence time is extremely short, we wondered whether tpe is a good species or not. If it is, populations of tpe are expected to be a monophyletic group. Two evolutionary scenarios were further tested following the founder speciation scenario: multiple origins (i.e., polyphyly) of tpe (Figure 2C) and a single origin followed by divergence (i.e., monophyly) in tpe (Figure 2D). The multiple origins scenario has a higher marginal density (5.577 × 10⁻⁴⁹, 99.836%) than the single origin (marginal density: 9.167 × 10⁻⁵², 0.164%), indicating a rejection of the monophyly of tpe. Similar to the estimation above, the effective population size of populations of tpe is smaller than that of bar (Figure 8E), while the divergence times have wide ranges from < 100 generations to 55,370 generations (95% CI: 6,814–861,759 generations; Figure 8F). The estimated mutation rate (μ) of each locus in multiple origins scenario has also been estimated (μc: 5.274 × 10⁻⁴, 95% CI: 1.999-8.9 × 10⁻⁴, Figure 8H). Such results of a wide range of divergence times of tpe populations from bar reflect (1) recent divergence and (2) sharing common ancestral polymorphisms that results in long coalescence times. This result also suggests that the divergence time is insufficient for complete speciation. Therefore, populations of both bar and tpe are combined into a single taxonomic group for the further analysis of local climatic adaptation. However, the very small migration rate on average (m < 0.001, Figure 8G) indicates that the gene flow seems to have ceased after population colonization.

After removing 13 bioclimative variables of collinearity, the remaining six bioclimatic variables are: bio2 [mean of monthly (max temp − min temp), named the Mean Diurnal Range], bio8 (mean temperature of wettest quarter), bio9 (mean temperature of driest quarter), bio13 (precipitation of wettest month), bio18 (precipitation of warmest quarter), bio19 (precipitation of coldest quarter). Among these six bioclimatic variables, the first three (bio2, bio8, and bio9) belong to the temperature dimension, while the last three (bio13, bio18, and bio19) are the precipitation dimension.

Although multicollinearity was prevented by removing variables with VIF > 10, certain variables were still partially intercorrelated. For example, the temperature variable bio2 was negatively correlated with the precipitation variable bio19, and two precipitation variables, bio13 and bio18, are negatively correlated (Supplementary Figure 3). The partial correlation of bioclimatic variables shown by PCA also revealed that the explanatory direction for the Grinnellian niche distribution is opposite for bio2 and bio19, and that bio13 and bio18 have similar explanatory intensity and directions (Figure 9). The PCA plot drawn on the first two axes explains 45.5 and 32.5% of the variation in the Grinnellian niche, but this biaxial PC distribution cannot separate bar and tpe well (Figure 9). When solely considering the temperature dimension (bio2, bio8, and bio9), or solely considering the precipitation dimension (bio13, bio18, and bio19), similar results of indistinguishable clustering between bar and tpe were revealed (Supplementary Figure 4). It is worth noting that the explanatory percentage of the first two PCs of the precipitation dimension (95.2%, Supplementary Figure 4B) was estimated to be higher than that of the temperature dimension (84.2%, Supplementary Figure 4A), suggesting that the precipitation could be more relevant in explaining the geographic distribution of these two species. Among these six non-collinear bioclimatic variables, only bio8 (mean temperature of wettest quarter) was shown by the MLGR analysis to be significant for predicting the occurrence of species (P = 0.036, Supplementary Table 2).

Ecological niche modeling (ENM) was further used to predict the suitable distributions of both bar and tpe as well as to examine the key climatic variables in the prediction. Distribution models of both species showed high discrimination performance. The cross-validation area under the curve (AUC) value for all models is 0.999, indicating that 99.9% of the records were correctly predicted. The climatic suitability of each species is proportional to similar climatic variables. Bio18 (precipitation of warmest quarter) and bio2 (mean diurnal range) are the first two contributors to the species distribution model in bar (57.9 and 28%, respectively) and in tpe (27.7 and 47.6%, respectively).

The ranges of both species predicted by the six bioclimatic variables (bio2, bio3, bio8, bio13, bio18, and bio19) are roughly consistent with the currently known distributions, except for a certain probability (< 50%) of suitable habitats predicted in the Hengchun Peninsula of southern Taiwan for the tpe (Supplementary Figure 1). In accordance with the specimen records, northern Taiwan is suggested as the most suitable habitat for both species by ENM. In contrast to the smaller ranges of tpe, the predicted distribution of bar is broader, including the most southern area and the northeastern and northwestern areas of Taiwan (Figure 1 and Supplementary Figure 1). The central area of Taiwan, where the Central Mountain Range is, is not suitable for either species distribution.

To assess whether the geographic or the environmental difference drives the genetic divergence among populations, isolation-by-distance and isolation-by-environment tests were conducted using the Mantel test. The Spearman correlation shows no significant correlation between geographic and environmental distances (ρ = 0.2888, P = 0.0872), between genetic and geographic distances (ρ = 0.2454, P = 0.0659), or between genetic and environmental distances (ρ = 0.1105, P = 0.3606; Table 4). The partial Mantel test did not detect significant correlation between genetic and environmental distance when conditioning on the geographic effect either (ρ = 0.2743, P = 0.0814, Table 4). Similar results of no autocorrelation between geographic and environmental distances (r² = 0.0152, β = 0.1201, P = 0.501) and non-significant effects of geographic (r² = 0.0373, β = 0.1418, P = 0.281) and environmental distances (r² = 0.0063, β = 0.0600, P = 0.685) on the change of genetic distance were obtained by MMRR analyses (Table 4). The joint effect of both geographic and environmental distances also did not affect the genetic distance significantly (r² = 0.0404, β_geo = 0.1367, P = 0.345, β_env = 0.0426, P = 0.797, Table 4). These results indicate that the population genetic differentiation (or gene flow) is not linearly correlated with the geographic and climatic differentiation.

Despite nonlinear correlation between the climatic differentiation and genetic distance, we still wondered whether the local climatic heterogeneity explains the genetic variation of populations. When conditioning on the geographic distribution, we found that 59.5% of the variation is explained by climatic variables (Table 5). The ANOVA further indicates that all six predictors (bio2, bio8, bio9, bio13, bio18, and bio19) can significantly explain population genetic components (P < 0.0001, Table 5), and the bio9 and bio2 have the highest explanatory proportions (23.0 and 11.7%, respectively) for predicting the population genetic variation (Table 5). The distribution of climatic variables along the ordination axis was further examined by GLM. Four bioclimatic variables including two variables (bio8 and bio9) of the temperature dimension and two variables (bio13 and bio18) of the precipitation dimension have significant F statistics (adjusted R² = 0.041, 0.104, 0.216, and 0.178, P < 0.05, Table 5), but only the precipitation variables bio13 and bio18 correlate significantly with the ordination axis1 of dbRDA, while all four variables are significantly correlated with axis2 (Table 5). However, the adjusted R² of bio8 could be too small (roughly explaining 4% of the variation on axis2) to be meaningful. If only the geographic clusters are considered (i.e., locations of populations) without genetic information, all the bioclimatic variables cannot predict the occurrence of populations by MLGR independently. Nevertheless, the joint effect of the four factors, bio8 × bio13 × bio18 × bio19, can predict the populations by MLGR significantly (χ² = 26.540, df = 8, P = 0.0008).

Scutellaria taipeiensis is a species morphologically alike to bar. It was named in 2003 based on its leaf shape (length-to-width ratio < 2 and triangular-ovate to broadly ovate shape) and the rounded concentric tuberculae type of nutlet coat that distinguished it from bar (Huang et al., 2003). However, in our SEM observation, we found that the nutlet coat pattern is not different between species and the “radiated umbrella-like” character (Hsieh and Huang, 1995) may disappear when seeds mature (Figure 3). In other words, the morphological differences of nutlet-coat patterns described by Hsieh and Huang (1995) and Huang et al. (2003) may result from the comparison of seeds at different developmental stages. In addition, the leaf size and shape is usually plastic in response to environmental changes (Rozendaal et al., 2006; Royer et al., 2009), and the range of leaf shape is wider and more varied in bar, covering the range of tpe (Supplementary Table 3). Morphological similarity in both vegetative and reproductive organs implies that tpe could merely be an ecotype or form of bar instead of a distinct species. To confirm this speculation, a comprehensive study of the genetic and ecological properties was carried out to complement the morphological evidence (Camargo et al., 2012).

The genetic variation of tpe is essentially a subset of that found in bar (Figure 6A), which matches the geographic distribution of these two species (Figure 1 and Supplementary Figure 1). In addition, both genetic and niche components are non-significantly different between species (Figures 7A–D, 9, Supplementary Figure 4). These analyses suggest that the speciation is incomplete, or at least that species divergence is not reflected in the overall genetic differentiation and niche preference. Our result resembles the common phenomenon that the progenitor species is paraphyletic to its derived species in plants because the ongoing gene flow increases the time for the derived species to achieve monophyly (Rieseberg and Brouillet, 1994). However, in this case, the daughter species tpe does not form a monophyletic group but fits the model of multiple founder speciation events (Figure 2B).

In fact, the estimated time of founder speciation of tpe is quite short (< 100 generations) in contrast to the long coalescent time of the Taiwanese populations of bar (11,145 generations ago, 95% CI: 1,505–92,795 generations ago, Figure 8D). If a generation is taken as 1 year, the coalescent time of the Taiwanese population of bar is roughly at the end of the Last Glacial Maximum (LGM), the time that Taiwan Island was most recently connected to continental Asia (Nino and Emery, 1961; Voris, 2000), which is probably the time when bar colonized Taiwan. However, the extremely short divergence time of tpe supports the hypothesis that tpe is not a good species but could be a subset of bar, or alternatively, that the loci that contribute to the population's differentiation are not involved in the species' divergence (Table 3). Shared polymorphisms across closely related species appear to be common in island species (Lopez-Sepulveda et al., 2013; Takayama et al., 2015), which might result from introgression following secondary contact (Minder and Widmer, 2008) or from maintaining ancestral polymorphisms due to incomplete lineage sorting within recently diverged species (Nagl et al., 1998). The ABC analysis suggested that sharing ancestral polymorphism could better explain the genetic structure of bar and tpe than secondary contact did. Furthermore, the best-fit scenario of multiple origins suggests insufficient time for lineage sorting in neutral genes after the divergence of speciation genes or speciation traits (i.e., the genes or traits involved in speciation). Multiple origins from a single progenitor usually are seen as an evolutionary mechanism for the recent and rapid radiation of island species rather than the mature radiation that sometimes leads to a high diversity of continental species (cf. Linder, 2008). However, such an evolutionary mechanism of multiple origins could be unstable to fusion and convergence to monophyly.

The Grinnellian niche that is comprised of temperature and precipitation is important and related to plant speciation, in particular in those budding speciation cases with peripheral ranges or distributional overlaps (Martinez-Cabrera et al., 2012; Grossenbacher et al., 2014). Although no positively selected loci was detected (Supplementary Figure 2), variation of neutral genes could still be affected due to demographic changes (e.g., population size decline, Figures 8A,B) under strong environmental pressures (Schoville et al., 2012). Bio8 is the only bioclimatic variable that has a significant effect on the prediction of groups of species. Bio8 is the mean temperature of the wettest quarter. However, a non-significant difference in bio8 was found between bar (24.72 ± 2.273°C) and tpe (24.66 ± 2.947°C) in means (independent two-group Mann–Whitney U-test, W = 132.5, P = 0.854) and variance (Kruskal–Wallis test, χ² = 0.039, df = 1, P = 0.844). This is probably because the overlapping distribution ranges of bar populations and tpe populations result in a similar climatic situation, thus causing a false positive for species prediction by the climate factor.

Despite the lack of species divergence in genetics and niches, local climatic heterogeneity has driven the genetic differentiation among populations within the island. This phenomenon is very like the case of Aeonium davidbramwelii and A. nobile on the island of La Palma in the Canarian archipelago: the former is an ecological generalist while the latter is an ecological specialist. There is a lack of species divergence but an obvious population structure exists (Harter et al., 2015). However, in our case, although the niche space of tpe is smaller than that of bar, their spatial distribution predicted by ENM on the basis of six bioclimatic variables is almost overlapping (Supplementary Figure 1). The only certain difference with the real distribution of tpe in the Taipei Basin is that the spatial model predicts a moderate probability for the potential distribution of tpe in southern Taiwan (Supplementary Figure 1). The dispersal mechanism of Scutellaria species is suggested to be attributable to rivers (Williams, 1992; Middleton, 2000; Chiang et al., 2012a). The lack of rivers running north to south in Taiwan decreases the probability of a southward long-distance dispersal of tpe.

Considering all genetic, evolutionary, ENM, and morphological evidence, we conclude that no species delimitation can be made between bar and tpe. Therefore, all populations are taken as one species. A climatic effect, in particular the precipitation dimension, leads to local population differentiation. Although, geographic and environmental distances do not matter in the genetic differentiation among populations (Table 4), we still propose that the climatic variables determine the genetic components of populations inferred by both partial dbRDA (Table 5) and MLGR significantly. These results suggest that climatic heterogeneity is not linearly correlated with the genetic composition, but could affect genetic structure independently and locally, i.e., climate leading to local adaptation.

The environmental heterogeneity, which causes significant genetic divergence (Figures 7E–H) and contributes a high proportion of population genetic variation (76.78%, Table 3), could be a major driver leading to postzygotic isolation (Martin and Willis, 2007; Anacker and Strauss, 2014). According to the AMOVA, roughly 3/4 of the variation can be attributed to population divergence (Table 3), and the partial dbRDA indicates that ~60% of the variation among populations can be explained by constrained variables (i.e., climatic factors, Table 5). This means that roughly half the genetic variation can be attributed to local climatic heterogeneity. Environmental heterogeneity could increase the genetic or genomic barriers among populations, and thus result in local adaption (Kawecki and Ebert, 2004; Savolainen et al., 2007).

In the analyses that did not take the genetics into account, the joint effect of bio8 × bio13 × bio18 × bio19 inferred by MLGR and of bio2/bio18 inferred by ENM are suggested as the main climatic variables affecting population clustering and distribution; if the genetic factor is considered, the bio13/bio18 on axis1 and bio9/bio13/bio18 on axis2 of the dbRDA are suggested as significant explanations for the genetic variation among populations. The climatic variables bio18 and bio13 are two key variables inferred from different statistical analyses (except bio13 by the ENM). These two variables belong to the precipitation dimension, standing for the precipitation of the warmest quarter and the precipitation of the wettest month, respectively (Supplementary Table 4). The significant effect of bio13 and bio18 was revealed in the multiple hills (clusters) of the ordisurf plots of bio13 and bio18 in partial dbRDA (Supplementary Figures 5D,E) in contrast to the other bioclimatic factors (Supplementary Figures 5A–C,F), and their high correlation is also revealed in the almost overlapping coordinate axes (i.e., arrows) of the original bioclimatic variables bio13 and bio18 (Supplementary Figure 5).

According to the WorldClim database, the wettest month is September (340.4 mm on average) and the warmest quarter of all recorded populations is July to September (29.63°C on average; Supplementary Figure 6). In fact, bio13 and bio18 are highly correlated (Spearman's rank correlation, ρ = 0.757, P = 1.75e–12). According to the common characteristics of bio13 and bio18 as well as their significant correlation, we suggest that the main determinant governing the genetic structure of populations is actually the precipitation from late summer to early autumn. This period is during the typhoon season (late summer and early autumn) in Taiwan and just after the flowering and fruiting season (tpe: Mar–Jun; bar: Oct–May; Hsieh and Huang, 1995; Huang et al., 2003) with subsequent germination and growth. This is the most important period for population regeneration in this group of species (i.e., bar and tpe). Precipitation during seed germination and growth could affect the demographic size and could also influence the success of seed colonization. In addition, since the seed dispersal of Scutellaria species is thought to be spread by water (Williams, 1992; Cruzan, 2001), the abundant rainfall in typhoon season may also contribute to the long-distance dispersal, and further affects the population genetic structure.

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.