2.1 Species and sampling locations

Mussaenda kwangtungensis H. L. Li (Rubiaceae) is a perennial, climbing shrub, between 1 m and 2.5 m in height. It has papery leaves with 3–9 cm × 1–3 cm in size and lanceolate-elliptic to elliptic-oblong in shape. Its inflorescences are compact cymose-to-subcapitate, and each possesses white, blade oblong-ovate, elliptic, or elliptic-ovate petaloid calycophylls with 3.5–5 cm long and 1.5–2.5 cm wide (Figure 1). The tubular flowers are yellow and 25–35 mm in length. This species flowers between May and September, with a peak between June and July; berry fruits mature between May and November (Luo et al., 2015). M. kwangtungensis is a typical functional dioecious species (Duan et al., 2018; Li, Wu, Zhang, & Barrett, 2010) which is pollinated by butterflies, moths, and bees (Luo et al., 2015). This species is endemic to China; it has a very narrow distribution, only occurring in South Guangdong Province (Figure 2) in secondary, disturbed forests, in open habitats, and on sunny slopes.

The study region assessed here encompasses the whole distribution range of M. kwangtungensis, comprising offshore islands in Guangdong Province, South China, as well as adjacent mainland areas (Figure 2). We collected samples from 23 populations (Table 1), including 12 island populations from 12 islands, alongside eight inland and three peninsula populations. The area of selected islands varies from 1.27 to 157 km² (Committee, 2013), and the nearest distance between an island and the mainland fell within a range of 3.8 to 77.7 km (Table S1 in Appendix S1). In each population, between eight and 33 individuals (mean = 27) were sampled, a total of 617 individuals. Each sample was separated from others by at least 30 m in order to reflect the real pattern of the population as much as possible. About five healthy leaves were collected from each individual and dried using silica gel. The voucher specimen for each population was deposited in South China Botanical Garden, Chinese Academy of Sciences (Table S2: Appendix S1).

2.2 DNA extraction, microsatellite genotyping, and chloroplast sequencing

Total genomic DNA of M. kwangtungensis was extracted from ca. 50 mg samples of silica gel-dried leaves using the modified cetyltrimethylammonium bromide (CTAB) method (Doyle, 1991). Thus, 11 nuclear microsatellite loci (i.e., CT99, CT142, CT48, CT59, CAA92, CT135, CT17, AC30, CT113, CT12, and CT60) were amplified for all 617 individuals using fluorescently labeled primers (Duan & Zhang, 2014; Appendix S1). Alleles were then scanned on ABI Prism 3100 Genetic Analyzer (Invitrogen, China) with internal size standard GeneScan™ 500 LIZ. and called by GeneMarker v. 2.4.0 (Holland & Parson, 2011). In order to minimize the genotyping error generated during different runs (Hoffman & Amos, 2005), the same laboratory methods, personnel, sequencing instrument, and allele bins in GeneMarker were used in all cases. Additionally, around 60 individuals were randomly selected as replicates, which were amplified and scanned at least twice in different runs.

For chloroplast sequencing, between two and eight individuals (mean = 5) in each population were successfully sequenced for three chloroplast intergenic spacer regions: ndhf-rpl32, psbE-petL, and rpl32-trnL (Shaw, Lickey, Schilling, & Small, 2007). Polymerase chain reactions (PCRs) were conducted following the protocols outlined in Appendix S1. PCR products were then sequenced by an ABI 3,730 automated DNA analyzer (Applied Biosystems). DNA sequences were assembled using Sequencher 5.3 (Gene Codes Corporation 2015), aligned with BioEdit 7.1.11 (Hall, 1999) and deposited in GenBank, under accession numbers MT349425–MT349442 and MT361334–MT361343.

2.3 Nuclear data analysis

2.4 CpDNA sequence analysis

2.5 Ecological niche modeling

To identify the potential species distribution during the LGM, we calibrated climatic envelope models using georeferenced native occurrence records for the species and the Maxent algorithm. Current presence records of M. kwangtungensis were collected from the herbarium records and our sampled populations (Table S2). Bioclimatic variables were downloaded from the WorldClim database (http://worldclim.org) for three different periods: the last interglacial (LIG, ~120–140 kya), the LGM (~21 kya), and current. Palaeoclimatic layers were drawn both from the Community Climate System Model (CCSM) (Collins et al., 2006) and the Model for Interdisciplinary Research on Climate (MIROC) 3.2 in simulations of the potential habitats during the LGM. A cross-validated model was then projected onto scenarios of LIG and LGM using MAXENT v.3.3.3 (Phillips, Anderson, & Schapire, 2006). In running models, we used the default values of convergence threshold (10^–5) and the maximum number of iterations (500), using repeatedly randomized samples of 25% of localities for cross-validation model training and testing, respectively. The regularization multiplier was set to 15 to fit smoother curves as these are thought to perform better when projected to different time periods, and because they do not overfit to training data. Model performance was evaluated using two different indices, the area under the curve (AUC) of the receiver operating characteristic and the true skill statistic (TSS). The predicted distribution at different periods was visualized on the map by ArcGIS v.10.2. Suitable areas for simulated models from MAXENT during different periods were calculated using ArcGIS by eligible considering distribution threshold greater than 0.5.

3.1 Nuclear genetic diversity and differentiation

Comparing with an adjusted P value, no linkage disequilibrium was found in any population among loci (p > .05) indicating that all 11 loci are independently inherited in M. kwangtungensis. Tests of deviation from the HWE, apart from CT12, revealed significant departures in four to 15 populations of other loci. Similarly, except for locus CT12, the other ten showed signals of null alleles in two to 15 populations with an average frequency between 3.15% and 13.83%. Analysis of genetic variation of M. kwangtungensis was therefore based on these 11 loci.

At the species level, totally 239 alleles from 11 microsatellite loci were amplified with a mean of 21.6 per locus. Overall genetic diversity (H_E) was 0.705. At the population level, the mean number of alleles (N_A) across all 11 loci ranged from 5.7 (BL) to 10.7 (SC). When correcting for the sample size, allelic richness (A_R) varied from 5.01 (WL) to 6.62 (SC). Observed heterozygosity (H_O) ranged from 0.494 (DX) to 0.680 (ZM) and expected heterozygosity (H_E) from 0.602 (DG) to 0.757 (BL) (Table 1). Of A_R, H_O, and H_E, none was significantly higher on the mainland than on islands (Table 3), indicating that genetic diversity was not significantly different between mainland and island populations. Correlating genetic diversity of island populations with the nearest distance to mainland showed no significant relationship (p > .08), while correlating with the area of islands, only allelic richness (A_R) was significantly positively correlated (p = .016) (Figure 3).

Overall F_ST across the 11 loci was 0.078 (p < .01), while standardized F'_ST was 0.247, indicating moderate differentiation. Genetic differentiation of mainland populations (F_ST = 0.073) was not significantly different compared to island populations (F_ST = 0.067) (p = .698) (Table 3). Our AMOVA revealed that only 0.90% of genetic variation resided among mainland and island populations (Table S4, Appendix S2). There was significant isolation by distance, as shown by the relationship of F'_ST/ (1 − F'_ST) with geographical distance (Mantel test) among all populations, as well as in both mainland and island populations (p < .05, Figure 4).

3.2 Nuclear genetic structure

Population structure and ancestry were inferred in STRUCTURE. The estimated logarithm of probability of the data, Ln P(K), increased slightly from K = 2, and the ΔK statistic was greatest for K = 2, with a second peak at K = 3 (Figure S2). At K = 2, all populations were divided into two clusters, geographically corresponding to inland and island plus peninsular regions, respectively. In terms of K = 3, inland cluster remained consistent, while DX, DG, XC, and SC were further clearly separated from island cluster, and WS, XW, and BL showed a large amount of mixture among gene pools (Figure 2).

3.3 Bottleneck and ABC analysis

Bottleneck analysis did not indicate excess of heterozygosity associated with recent bottlenecks in any mainland and island population, either under the SMM or TPM models (Table S5: Appendix S2). However, a shift from the normal L-shaped distribution of allelic frequencies was found in the BL case.

A comparison of posterior probabilities of the seven scenarios using a logistic approach on 70,000 data points indicated that scenario 2 had the highest support (Figure S3), with a probability of 0.457 (95% CI: 0.365–0.549), associated with the hypothesis that there was one refugium in island and that mainland populations were derived from this area. Scenario 1 had the second highest support with a probability of 0.337 (95% CI: 0.229–0.447). Posterior probabilities for scenarios 3, 4, 5, 6, and 7 were 0.060, 0.053 0.023, 0.061, and 0.007, respectively. Posterior distribution of parameters for scenario 2 is shown in Table 4. Thus, if we consider a mean generation time of 10 years for M. kwangtungensis, the divergence of inland cluster (gene pool 1) from island cluster (gene pool 2), designated by t1, appears to have started 3,930 BP (95% CI: 1,020–7,140) compatible with the time of island formation during the Holocene. The median value of t2 (the time when ancestral effective population size varied) was 61,200 BP (95% CI: 8,120–96,700), indicating an increase of effective population size during the LGM (Table 4).

3.4 Contemporary and historical gene flow

Based on BAYESASS analysis, the contemporary migration rate (m) from the mainland to islands was 0.005 (−0.002 to 0.013) while that from islands to the mainland was 0.016 (0.004–0.028). The overlapping 95% confidence intervals indicated no asymmetrical gene flow between island and mainland clusters. Similarly, MIGRATE did not reveal a pattern of asymmetrical gene flow between island and mainland. The relative migration rates (M) computed by MIGRATE were 4.450 (3.767–5.067) from the mainland to islands and 3.317 (2.533–4.033) in the other direction. When converting M to m for comparison, the median value of historical migration rate from the mainland to islands was 0.000445–0.00445 with μ ranging from 10^–4 to 10^–3, nonsignificantly different from contemporary migration rare, and vice versa.

3.5 CpDNA structure

By aligning, M. kwangtungensis sequences yielded lengths of 920 bp, 1,224 bp, and 1,386 bp for rpl32-trnL, psbE-petL, and ndhf-rpl32, respectively. A total of 28 variable sites, comprising 26 point mutations and two indels, were detected, defining 21 haplotypes (Table S6: Appendix S2). The geographical distribution of haplotypes is shown in Figure 5. The most widely distributed haplotype (H1) was found in 12 populations (52.2%) and 36 individuals, followed by H7 occurring in seven populations and 20 individuals. Six haplotypes (H1–3, H7, H8, and H12; 28.6%) were shared among mainland and island populations, while six (H4–6 and H9–11, 28.6%) and nine (H13–21, 42.9%) were privately restricted to mainland and islands, respectively. The parsimony network revealed closely related haplotypes (Figure 5b). The shared ones (H1, H7, and H3) were the most interior nodes and could be considered as ancestral. All the unique haplotypes in mainland and islands differed from them by one to three mutation steps, except for H11 by six steps from H3.

For M. kwangtungensis, within-population haplotype diversity (h) was 0.471, and nucleotide diversity (π) was 0.370 × 10^–3. On the mainland, each population had 1 (BZ) to 4 (ZM) haplotypes with a mean = 2.5. The value of h varied from 0 to 0.867 with a mean of 0.558, while π ranged from 0 to 1.327 × 10^–3 (YM) (Table 1). On islands, the number of haplotypes per population varied from 1 (SC, BL, DG, WL, SM and XX) to 5 (XW), while values for h ranged from 0 to 1.000 (XC and WS) with a mean of 0.391 and π ranged from 0 to 1.137 × 10^–3 (XC). None of diversity parameters (h and π) showed significant difference between mainland and island populations (p > .269). Correlation analysis between diversity at cpDNA (h and π) and island characteristics (area and distance to mainland) revealed no significant relationship (p > .160).

Population differentiation in cpDNA was revealed by F_ST = 0.402. Island populations showed a slightly stronger genetic differentiation (F_ST = 0.459) than their mainland counterparts (F_ST = 0.315). The results of AMOVA showed that just 0.34% of genetic variation existed among mainland and island groups (Table S4). No pattern of isolation by distance was revealed by Mantel tests in either all populations or within each group (p > .201).

3.6 Demographic history

At the species level, the neutrality test showed that both Tajima's D (−1.696, p = .012) and Fu's Fs (9.885, p = .003) were negative. The mismatch distribution for haplotypes was unimodal with nonsignificant SSD and H_Rag (p > .05, Table 5), which was statistically consistent with the expansion model. Based on the corresponding τ value, the demographic expansion event was estimated to have occurred at 21,200 (10,700–32,200) years BP. Additionally, the BSPs analysis detected a slightly increasing population size until 20,000 years BP, then following a recent decrease (Fig. S4). When analyzed separately, island cluster showed the same pattern: negative D (−1.420, p = .052 and Fs (−5.133, p = .022), unimodality of mismatch distribution (Table 5), and a slightly ascending curve followed by a fall for BSP (Fig. S4). The expansion event was dated at 17,000 (8,640 – 25,900) years BP. No signature of demographic expansion was detected in the case of the inland cluster. Both D (1.032, p = .810) and Fs (1.723, p = .837) were positive. The mismatch distribution was bimodal with significant SSD (Table 5). BSP showed a very slightly ascending trend (Figure S4).

3.7 Ecological niche modeling

The prediction of the bioclimatically suitable areas for M. kwangtungensis during the LIG, the LGM, and current is shown in Figure 6. The models in each period were indicated to have high predictive abilities with AUC (>0.99) and TSS (0.71–0.85). It is clear that M. kwangtungensis has been narrowly distributed at all stages. In general, the predicted distribution under current conditions was similar to the actual distributions. CCSM model suggested the habitats available during the LGM were relatively limited compared with the MIROC model. However, both of them showed the same trend: M. kwangtungensis experienced a recently spatial expansion when temperature decreased during the LGM period. In particular, the area with climatically suitability >0.5 enlarged from the present distribution to the southern ancient coast. Thus, by calculating the area with climatic suitability >0.5, larger region was highly suitable for M. kwangtungensis during the LGM based on CCSM (17,183 km²) and MIROC (29,560 km²) models, compared with that at present (10,118 km²). The area with climatic suitability >0.5 from the LIG (18,433 km²) was a little larger than that during the LGM based on the CCSM model, but far smaller than that from LGM based on MIROC model. Data show that when climatic suitability was greater than 0.8, CCSM model revealed greater area in inland than on islands, while in MIROC model, more suitable habitats were located on islands.

4.1 Similar levels of genetic diversity between mainland and island populations

Plastid and nuclear markers reveal high levels of genetic diversity in this species. Indeed, haplotype diversity based on cpDNA (h_T = 0.855) was higher than the mean value of angiosperms (h_T = 0.68) (Petit et al., 2005), while the expected heterozygosity computed on microsatellites (H_E = 0.705) was higher than values reported for narrow-distributed/endemic species (H_E = 0.56/0.42) (Nybom, 2004). Life-history traits, especially mating systems, may determine genetic structure among populations (Duminil et al., 2007). Species with mixed or outcrossing mating systems was demonstrated to have significantly higher genetic diversity than selfing species (Nybom, 2004). It is clear that M. kwangtungensis is a typical functional dioecious species with an obligate outcrossing mating system (Li et al., 2010), thus considered capable of harboring high genetic diversity (Hamrick & Godt, 1996). This has been evidenced by many previous studies (Shi & Chen, 2012; Shi et al., 2014).

Island populations show comparable genetic variation with mainland populations (Table 3), inconsistent with the general expectation that these tend to have lower genetic variation due to founder effects and geographical isolation (Frankham, 1997). We discuss possible explanations for this as similar, or higher, levels of genetic variation on islands have been variably ascribed to a number of variables including high gene flow. In this first case, we know that M. kwangtungensis is pollinated predominantly by butterflies (Papilionidae: Papilio paris, Figure 1) (Luo et al., 2015), in which long-distance dispersals have been frequently observed (up to 13 km) (Baguette, 2003; Li, Zhang, Settele, Franzén, & Schweiger, 2013). Individuals of P. paris have been estimated to move up to 6 km (personal communication with Dr. X. Li), a distance that ensures pollen can be transferred between most islands and adjacent populations (Table S1). In addition, the color of berry fruits turns to purple black when they mature and attracts foraging birds (personal observation), which may provide opportunities for the long-distance dispersal of seeds. At the same time, however, the high levels of private haplotypes on islands (Figure 5) also seem to indicate a lower level of occurrence while similarity in contemporary and historical gene flow indicate that isolation does not exert an influence on genetic structure (Chiucchi & Gibbs, 2010). Considering a relatively continuous distribution on the mainland, sufficient dispersal across the ocean barrier could be expected. Secondly, genetic bottlenecks on the mainland might provide another explanation as these have occurred in the recent past due to human disturbance. The construction of buildings and roads, for example, has been reported in mainland populations, including for Ficus pumila (Chen et al., 2008). Bottleneck analysis here, however, did not show any evidence for a population bottleneck on the mainland, as well as on islands excepting for BL (Table S5), indicating a relatively stable effective population size, as revealed by neutrality tests and mismatch analysis based on cpDNA. Thirdly, multiple mainland–island colonization may have occurred given that species originated on the mainland and subsequently colonized islands. This would have led to a decrease in genetic diversity with increased distance to the mainland due to the founder effect under a stepping-stone colonization process (Yamada & Maki, 2012). However, our results did not indicate any correlation between any parameter of genetic diversity and distance to mainland. The haplotype network also did not reveal a clear divergence between mainland and islands; ancient haplotypes (H1, H3 and H7) were distributed in both regions (Figure 5) and this does not support mainland–island colonization. Fourthly, consideration of islands as refugia during the LGM (Chen et al., 2008; Fernández-Mazuecos & Vargas, 2011; Tsumura & Ohba, 1993) can be based on relict populations frequently reported to preserve high genetic variations (García-Verdugo et al., 2013; Hampe & Petit, 2005; Petit et al., 2003). In one example, Cryptomeria japonica maintained a relatively large genetic variation on Yaku Island compared with that of mainland Japan mainly because this island was refugium during the LGM and had a suitable climate (Tsumura & Ohba, 1993). In the case of the polyploid species Prunus lusitanica, island populations displayed more private alleles and higher genotypic diversity in old volcanic areas, suggesting the idea of oceanic islands as important refugia for biodiversity (García-Verdugo et al., 2013). Our data also supported the hypothesis of islands as refugia in a very young, continental island system (detailed discussion seen below).

4.2 No effect of isolation on genetic diversity and differentiation

As proposed by the equilibrium theory of island biogeography (Losos, Ricklefs, & MacArthur, 2009; MacArthur & Wilson, 1967), larger islands can be expected to hold higher cumulative levels of diversity because of larger numbers and sizes of populations. This means that levels of genetic diversity on islands should correlate with island area (McGlaughlin et al., 2014). However, we found that island area was only a predictor for one measure of genetic diversity (A_R) based on microsatellites in M. kwangtungensis (Figure 3), and had no significant effect on other measures either for microsatellites or cpDNA. These findings are similar to those found for two endemic island species of Acmispon on California Channel Islands (McGlaughlin et al., 2014), indicating that the correlation between genetic diversity and island area is not a common pattern. One possible explanation is that although a larger island is capable of sustaining more populations of a species, area should not be taken as an accurate surrogate for effective population size as this will be more significant determinant of levels of genetic diversity than number of populations (García-Verdugo et al., 2013; McGlaughlin et al., 2014). The prediction of a positive relationship between genetic diversity and island area assumes that population densities for a species are kept the same across different islands, something that is actually unlikely to be the case. Indeed, for M. kwangtungensis, XX island with the smallest area was observed to have an obviously higher population density, compared to other larger islands (personal observation).

A negative relationship between distance to mainland and genetic diversity is predicted, mainly due to founder effects and limited gene flow. Although this was also seen in two previous studies on Izu Islands of Japan (Inoue & Kawahara, 1990; Yamada & Maki, 2012), we failed to establish a significant relationship in the case of M. kwangtungensis. The most likely explanation for this result is that island systems had different histories of formation and species colonization. In the case of the Izu Islands, these are distributed as a roughly linear chain extending from north-to-south near Mainland Japan and are thought to be formed southwards during the late Pleistocene through submarine volcanism (Kaneoka, Ishiki, & Zashu, 1970; Yamada & Maki, 2012), probably leading to a southward stepping stone-like colonization from the mainland for Weigela coraeensis (Yamada & Maki, 2012) and Campanula punctate (Inoue & Kawahara, 1990). The islands studied here, in contrast, are distributed along the coast of the South China Mainland, which were connected during the LGM and were formed by rising sea levels around 7,000 years ago (Chen et al., 1993). Indeed, M. kwangtungensis is thought to have survived long term in these islands based on genetic and ecological niche modeling (Figure 6, see discussion below), rather than as a colonization/recolonization from a mainland source after the ice age.

Restricted land masses of islands with the geographical barrier of the sea are inclined to promote genetic differentiation among islands, and between the mainland and island populations (Frankham, 1997). Our results revealed statistically nonsignificant differences between genetic differentiations of mainland and islands at microsatellites (Figure 4, Table 3) as well as low levels of genetic differentiation between them, as indicated by the results of AMOVA for both genetic markers (Table S4). These results may be attributable to potentially high rates of gene flow (see discussion above) as well as the geographical proximity of these areas. The islands assessed here are closely distributed along the coastal line with distances to mainland ranging from 3.8 to 77.7 km and with nearest neighbors from 0.7 to 19.6 km. Thus, SM, DJ, and XX islands are within a range <6 km from the nearest peninsula (Table S1); this close proximity is likely to give rise to a high rate of gene flow, as suggested by patterns of isolation by distance (Figure 4), and the STRUCTURE analyses based on microsatellite data, where these three islands with peninsular populations (QN, YM, and DX) fell into the same cluster (Figure 2).

4.3 In situ glacial survival on islands

Consistent with simulated vegetation maps (Harrison et al., 2001) and sporopollen evidence (Chen et al., 1993) that a land bridge in South China provided refugia for plant species during the LGM, our findings also suggest in situ glacial survival of M. kwangtungensis on islands. In the first place, ENM indicates a large area with highly climatic suitability was present in regions of this land bridge during the LGM (Figure 6). Although CCSM predicted relatively limited region compared with MIROC, both verified highly suitable habitats (>0.8) on islands, implying that the current populations are partly representative of the past. Secondly, molecular evidence revealed high genetic diversity and, especially, numbers of unique haplotypes within island populations (Figure 5, Table 1). In one example, BL showed the highest level of expected heterozygosity, with a minimum sample size of eight, while XW possessed the most haplotypes and, XC and WS had the maximum haplotype diversity (Table 1). Thirdly, STRUCTURE analysis based on microsatellite data revealed two clearly genetic clusters largely separating island populations from their mainland counterparts (Figure 2), indicating that islands represent a differentiated gene pool. Based on the clusters identified by STRUCTURE, our ABC analysis indicated the highest support for the hypothesis of one refugium in islands (scenario 2) (Table 2, Figure S3), although scenario 1 cannot be completely rejected due to overlapping confidence intervals. Fourth, network of cpDNA haplotypes revealed that H7, widely distributed on islands and the mainland, were located into the interior node (Figure 5), representing the ancestral one. These patterns collectively suggest islands have acted as important refugia for M. kwangtungensis during the LGM and that extant island populations are the products of in situ glacial survival.

It is the case that in situ glacial refugia have been frequently reported for temperate (Parisod & Besnard, 2007; Zeng et al., 2015) and evergreen broad-leaved forests (Chen et al., 2020; Gong et al., 2016; Wang et al., 2015), and mountainous regions have usually been regarded as critical glacial refugia because their high topographical heterogeneity offers scope for shifts in altitudinal range in response to climate changes (Shi et al., 2014). Islands, due to their stable climate and environmental heterogeneities, also have the potentials for in situ glacial survival of species, as evidenced by accumulating studies (García-Verdugo et al., 2013, 2019; Hutsemékers et al., 2011). In South China, land bridge, connecting between continental islands and mainland during glacial periods, may provide suitable habitats with a minor shift of temperature and humidity (Chen et al., 2012; Weigelt et al., 2013), and have supported effective population sizes for in situ persistence, as indicated in Hainan Island (Chang et al., 2012) and Taiwan Island (Cheng, Hwang, & Lin, 2005).

4.4 Demographic history

Typical responses of plants to climate changes include geographical migration, resulting in the alteration of distribution range (Etterson & Shaw, 2001). The current distribution of M. kwangtungensis is mainly restricted within a narrow belt in southern Guangdong Province. Based on ENM reconstruction, the possible occurrence of populations during the LGM mainly covered a greater area extending current distribution to southern land bridges (Figure 6), indicating much more suitable habitats for survival. Thereby, the growth of population size is considered to be expected. Our genetic evidence supported this idea: significantly negative D and Fs, unimodality of mismatch distribution (Table 5), and a slightly ascending curve for BSP (Figure S4) for all populations and island cluster were evident for population expansion (10,700–32,200 years BP), while inland populations kept a relatively stable population size. It is plausible that when climate cooling, land bridge connected islands to mainland, facilitating island populations expanding in space and population size, which may appear to be an unusual scenario, different from the general pattern of demographic contraction in glacial periods, but population growth in response to glaciation-induced increase in lowland area has recently been suggested for wide variety of species (Chen et al., 2012; Nakamura et al., 2014).

ABC analysis indicated the highest support for the hypothesis that there was one refugium in islands and inland populations were derived from island at around 3,930 years ago, when considering a generation time of 10 years for M. kwangtungensis (scenario 2) (Table 2, Figure S3). This scenario should be treated with caution given the overlapping confidence intervals with scenario 1 (one refugia in inland and island derived from inland). Since hypothesis of more than two refugia was rejected (Table 2, Fig. S3), one continuous refugia, where inland and island situated, may be an acceptable hypothesis, as predicted by both CCSM and MIROC during LGM (Figure 6). Both scenario 1 and scenario 2 supported a divergence time of ca. 3,000–4,000 years ago, compatible with the time of island formation, implying that vicariance sundered to the continuous distribution range of ancestral lineages (Nakamura et al., 2010).