Multilocus phylogeography and ecological niche modeling suggest speciation with gene flow between the two Bamboo Partridges

Abstract

Background

Understanding how species diversify is a long-standing question in biology. The allopatric speciation model is a classic hypothesis to explain the speciation process. This model supposes that there is no gene flow during the divergence process of geographically isolated populations. On the contrary, the speciation with gene flow model supposes that gene flow does occur during the speciation process. Whether allopatric species have gene flow during the speciation process is still an open question.

Methods

We used the genetic information from 31 loci of 24 Chinese Bamboo Partridges (Bambusicola thoracicus) and 23 Taiwan Bamboo Partridges (B. sonorivox) to infer the gene flow model of the two species, using the approximate Bayesian computation (ABC) model. The ecological niche model was used to infer the paleo-distribution during the glacial period. We also tested whether the two species had a conserved ecological niche by means of a background similarity test.

Results

The genetic data suggested that the post-divergence gene flow between the two species was terminated before the mid-Pleistocene. Furthermore, our ecological niche modeling suggested that their ecological niches were highly conserved, and that they shared an overlapping potential distribution range in the last glacial maximum.

Conclusions

The allopatric speciation model cannot explain the speciation process of the two Bamboo Partridges. The results of this study supported a scenario in which speciation with gene flow occurring between the allopatric species and have contributed to our understanding of the speciation process.

Background

How species multiply has been a fundamental question in evolutionary biology since Darwin published “On the Origin of Species” (Darwin 1859; Mayr 1942; van Doorn et al. 2009; Ronco et al. 2021). Understanding the speciation process is important in evolutionary biology to help us identify the genetic basis of local adaptation (Kautt et al. 2020), find the genetic loci linked with phenotype traits (Lamichhaney et al. 2015), and understand the biodiversity form process (Wang et al. 2018). Allopatric speciation supposes that geographically separated populations have no gene flow between them due to a geographical barrier, and that the populations are reproductively isolated due to genetic drift or other evolutionary forces (Mayr 1942, 1963; Wright 1943). Several lines of evidence support the allopatric speciation model (Mayr 1954; Lessios 1984; Coyne and Orr 2004; Boucher et al. 2016); however, there is also evidence for the alternative model, speciation with gene flow (Nosil 2008; Kautt et al. 2020; Wang et al. 2020).

The speciation with gene flow model holds that a geographical barrier is not a necessary condition for reproductive isolation, but that the speciation process can be completed even with gene flow taking place (Pinho and Hey 2010). Different ecological niches or gradual environmental changes can reduce gene flow between populations in different locations until, finally, there is reproductive isolation. The speciation with gene flow model is supported by some evidences, with examples including the divergence process of sympatric hawthorn and apple host races of Rhagoletis pomonella (Bush 1969; Egan et al. 2015), sympatric neotropical cichlid fishes (Amphilophus spp.) (Kautt et al. 2020), and parapatric shorebird species (Charadrius alexandrinus and C. dealbatus) (Wang et al. 2019). Sympatric species and parapatric species may not have had significantly different ecological niches in the initial stages of speciation or there may also have been dispersal between young species. The question, however, is whether allopatric species have gene flow between them during the speciation process.

Allopatric species of organisms occurring on continental islands and close neighboring species on the adjacent mainland may have no gene flow with each other during the divergence process (Coyne and Orr 2004). However, the continental islands would have had a “land bridge” during the glacial period (Voris 2000), which may have provided an opportunity for ancestral populations to share the same habitats and facilitate gene flow with each other during secondary contact. Fluctuations in sea level can cause alternation between allopatric and sympatric distribution, as seen throughout evolutionary histories. It is therefore difficult to distinguish allopatric speciation from speciation with gene flow for allopatric species, especially when the allopatric species are isolated by a strait. Therefore, we explored whether the allopatric species distributed on the mainland and continental islands showed evidence of gene flow during the speciation process.

The allopatric species pair, the Chinese Bamboo Partridge (Bambusicola thoracicus) and the Taiwan Bamboo Partridge (B. sonorivox), which is geographically isolated by the Taiwan Strait, could serve as an ideal model for testing whether there is evidence of gene flow during the divergence process. Both B. thoracicus and B. sonorivox are small partridges commonly occurring at the edge of evergreen broadleaf forests in south-eastern mainland China and Taiwan Island, respectively (Hung et al. 2014; Zheng 2015). The plumage coloration of the two Bamboo Partridges is highly distinctive (McGowan and Madge 2010), while their vocalizations are similar but distinguishable between species (Hung et al. 2014). A recent genetic study suggested that the two species diverged 1.8 million years ago based on the isolation with gene flow (IM) model (Hung et al. 2014). However, the IM model assumes that post-divergence gene flow was present during the entire divergence history of the two Bamboo Partridges (Strasburg and Rieseberg 2010). It fails to provide information about how long post-divergence gene flow has persisted and the corresponding level of gene flow (Strasburg and Rieseberg 2010; Hung et al. 2014). The lack of such information greatly hinders our understanding of the process and mechanism of their speciation.

In the present study, we combined genetic data, ecological niche modeling, and morphometrics to explore the speciation demography and level of ecological divergence between B. thoracicus and B. sonorivox, with the aim of inferring the speciation process. First, we inferred the duration and magnitude of post-divergence gene flow between the two Bamboo Partridge species. Using ecological niche models (ENMs), the fundamental ecological niches (the set of abiotic variables in which a species can persist, the variables include temperature and precipitation) were then estimated based on presence data of the species. With these analyses, we aimed to find evidence of gene flow during the speciation process of allopatric species, thus deepening our understanding of the speciation process.

Methods

Tissue sampling and DNA extraction

We collected blood and tissue samples from 24 individuals of B. thoracicus across their range in mainland China, and 23 individuals of B. sonorivox from Taiwan Island, China (Fig. 1; Additional file 1: Table S1). Samples were preserved in 75% alcohol and stored at − 80 °C for long-term storage. We used a DNA extraction kit (TianGen Biotech, Beijing, China) to extract gross genomic DNA following the manufacturer’ s protocol. DNA samples were stored at − 20 °C until later use.

Fig. 1

Sample distribution. Squares indicate the location of B. thoracicus sampling. Circles mark the location of B. sonorivox sampling.

Nuclear locus amplification, sequencing, and polymorphism

A total of 31 unlinked autosomal exon loci (Liu et al. 2018) and mitochondrial cytochrome b genes were amplified by polymerase chain reaction (PCR). The details of the PCR protocol and primer information are listed in Additional file 1: Table S2. Amplicons were sequenced on an Applied Biosystems 3730 capillary sequencer. Sequences were manually checked with the aid of DNAMAN 8.0.8 (http://www.lynnon.com) and aligned using Sequencher 4.7 (Gene Codes Corporation, MI, USA). We used the software PHASE 2.1 (Stephens et al. 2001; Stephens and Scheet 2005) which implemented in DnaSP 5 (Librado and Rozas 2009) to resolve diploid genotype sequences into haplotypes. Only haplotypes with PHASE probabilities higher than 60% were used for subsequent analyses. We calculated 7 intra- and interspecific DNA summary statistics for polymorphisms of each nuclear gene (Additional file 1: Table S3). The results were used for the subsequent approximate Bayesian computation analysis (ABC) (Beaumont et al. 2002). We followed the published method (Li et al. 2010) to estimate the substitution rate of each nuclear gene, and set the generation time of the Bamboo Partridge as two years.

Inference of speciation process for two Bamboo Partridges

We conducted an ABC analysis to fit the observed polymorphic pattern in the two Bamboo Partridges. To elucidate the details of the diversification history of the two Bamboo Partridges, we proposed two scenarios of post-divergence gene flow (Fig. 2). First, due to the small amount of gene flow estimated by the IM model (Hung et al. 2014), we assumed that post-divergence gene flow was only present during the early stage (the first 20%) of their divergence history. Alternatively, we proposed that the post-divergence gene flow could have lasted much longer (greater than 20% of their divergence history). For both models, the termination of post-divergence gene flow could not have occurred more than 20,000 years ago as the sea level between Taiwan and the Asian continent rose, thus interrupting any potential gene flow between the ancestor of the B. thoracicus and B. sonorivox. In contrast to the divergence with the gene flow model, we also used the isolation model as the null model, assuming no gene flow between the populations since their divergence. In all models, we estimated the population divergence time (T₁), termination time of post-divergence gene flow (T₂), long-term effective population sizes for both Bamboo Partridges (N_Bt and N_Bs for B. thoracicus and B. sonorivox, respectively), and the effective population size of their most recent common ancestor (N_a). For the two post-divergence gene flow models, we estimated the number of immigrants per generation: M_Bt and M_Bs, for the B. thoracicus and B. sonorivox, respectively, (M_Bt = 4 N_Bsm_Bt and M_Bs = 4 N_Btm_Bs, where m_Bt and m_Bs indicate the proportions of immigrants within B. thoracicus and B. sonorivox in each generation).

Fig. 2

Three demographic scenarios simulated by msABC. a Isolation model; b Early gene flow model; c Late gene flow model. N_a: Effective population size of the most recent common ancestor. N_Bt: Effective population size of B. thoracicus. N_Bs: Effective population size of B. sonorivox. m_Bt: Migration rate from B. sonorivox to B. thoracicus. m_Bs: Migration rate from B. thoracicus to B. sonorivox. T₁: Divergence time between B. thoracicus and B. sonorivox. T₂: Time of gene flow disappearance. PP: Posterior probability of each model.

We used the program msABC (Pavlidis et al. 2010) to perform coalescent simulations for all neutral loci. Demographic parameters were sampled from a uniform random distribution within their prior ranges. The ranges of the priors for demographic parameters in all models and the command line for each model are listed in Additional file 1: Tables S4 and S5, respectively. One million simulations were performed for each model, assuming the same sample size as the empirical data for each locus. Because the lengths of the loci used in the current study are short (< 1.5 kb), and linkage disequilibrium usually decays slowly in birds (Kawakami et al. 2017), we assumed that there was no intra-locus recombination in coalescent simulations. Furthermore, we used single-nucleotide polymorphism (SNP) frequency-based summary statistics in our ABC analysis that are relatively robust for demographic inferences when the lack of intra-locus recombination assumption is violated (Ramírez-Soriano et al. 2008). We scaled all time estimates to the unit of four N_Bt generations. The seven summary statistics (T₁, T₂, N_Bt, N_Bs, N_a, M_Bt, and M_Bs) mentioned above were calculated for coalescent simulation.

The program msABC (Pavlidis et al. 2010) was used to perform model selection and to infer the population parameters of the preferred model for the ABC analysis. We set the tolerance of the ABC analysis to 0.05. A multinomial logistic regression method (Beaumont et al. 2002) was used to infer the posterior probability of each demographic model for model selection. Then, the results of four million additional simulations were used to estimate the population parameters of the selected models using the neural network method (Blum and François 2010). We set the number of neural networks to 50. The demographic parameters were logit transformed, except for T₁, which was log transformed. This ensured that the posterior densities fell within the ranges of the prior distributions. We used the Epanechnikov kernel to calculate the posterior densities of each demographic parameter. We reported the median, mean, and mode of each parameter and the range of its 95% highest probability distribution (HPD) interval. The medians were used as point estimates of the demographic parameters.

Ecological niche models of the two Bamboo Partridges

To obtain the proxy of fundamental niches for both B. thoracicus and B. sonorivox, we generated their ENMs using Maxent v3.3.3 (Phillips and Dudík 2008). A total of 364 occurrence data points for B. thoracicus (Additional file 1: Table S6) were initially obtained from the Global Biodiversity Information Facility (GBIF) dataset (http://www.gbif.org) and birdwatching records (http://www.birdreport.cn/). For B. sonorivox, we retrieved 1162 occurrence data points (Additional file 1: Table S6) from the eBird dataset (https://ebird.org/). The study area included the minimum convex polygon surrounding the original occurrence records for B. thoracicus with a buffer of four decimal degrees; however, Taiwan and Hainan Island were excluded because of the absence of B. thoracicus on these large continental islands and the limited dispersal ability of these birds (Additional file 2: Figure S1). For B. sonorivox, only the entirety of Taiwan Island with a one decimal degree buffer was included in the ENM analysis (Additional file 2: Figure S2). The spatial resolution of the ENMs was 30 arc-seconds, c. 0.9 km, for both species. To prevent ENMs being biased by over-sampling, we randomly chose points that were at least 10 km and 1 km apart from one another from the original data points of B. thoracicus and B. sonorivox, respectively, which resulted in 351 occurrences for B. thoracicus (Additional file 1: Table S7; Additional file 2: Figure S1) and 568 for B. sonorivox (Additional file 1: Table S7; Additional file 2: Figure S2). We used 9 of the 19 bioclimatic variables from WorldClim (Hijmans et al. 2005; Fick and Hijmans 2017) as the environmental data layers for the ENMs by removing multicollinearity using Pearson’s correlations with r >|0.9| (Warren et al. 2010) (Additional file 1: Table S8; Additional file 2: Figure S3). The selected variables were mean diurnal range (Bio 2), isothermality (Bio 3), temperature annual range (Bio 7), mean temperature of the driest quarter (Bio 9), mean temperature of the warmest quarter (Bio 10), annual precipitation (Bio 12), precipitation of driest month (Bio 14), precipitation seasonality (Bio 15), and precipitation of the warmest quarter (Bio 18). We randomly set 75% of the data for training; the remaining 25% of the data were used for testing. To construct the ecological niche model with Maxent, we set the replicated run type as subsample, maximum iterations as 5000, apply threshold rule as minimum, all samples to background, and replicates as 100. The area under the receiver operating characteristic curve (AUC) value was used to evaluate the performance of the ENMs. For both Bamboo Partridges, if the suitability score of each cell generated by Maxent was larger than the mean minimum training presence logistic, it would be classified as suitable; otherwise, it would be considered unsuitable. We used two methods implemented in Maxent, namely percent contribution and permutation importance, to evaluate the relative contribution of each climatic variable to the ENMs of both Bamboo Partridge species.

Ecological niche comparisons

We quantified the differentiation (or overlap) between ecological niches of B. thoracicus and B. sonorivox based on their ENMs with ENMTools (Warren et al. 2010) using Schoener’s D (Schoener 1968) and similarity index I (Warren et al. 2008). Both indices were close with one when the ecological niche were more similar, in contrast, both indices were close with zero when the ecological niche were more different. We conducted a background similarity test (Warren et al. 2008) to test whether the actual ecological niche was more similar compared with the stochastically simulated ecological niche. For this test, we compared the empirical Schoener’s D and similarity index I values to a distribution of 200 stochastically simulated values. The stochastically simulated values were generated using the ENMs of one Bamboo Partridge and an ENMs created with random points drawn from the minimum convex polygon surrounding the original occurrence records of the other Bamboo Partridge. The geographic ranges of B. thoracicus and B. sonorivox were used as backgrounds individually. For the background test, we used the minimum convex polygon surrounding the original occurrence records and the entire Taiwan Island as the species range of B. thoracicus and B. sonorivox, respectively.

Potential suitable habitable regions of the two Bamboo Partridges during the last glacial maximum

To infer the likelihood of historical secondary contact between the B. thoracicus and B. sonorivox during previous glacial periods when the sea level was low, we reconstructed their potential distribution ranges in the late Pleistocene by projecting their ENMs on most parts of East Asia (Additional file 2: Figure S3) using simulated palaeo-climate data for the last glacial maximum (c. 22,000 years ago). We used bioclimate variables generated by the GCM model at 2.5 arc-minutes, c. 4 km, resolution for the last glacial maximum paleoclimate.

Estimation of the differentiation of morphology

We measured the morphological characteristics of 24 specimens of B. sonorivox (10 males and 14 females) and 23 of B. thoracicus (16 males and 7 females) archived in museums and research institutes (Additional file 1: Table S9). Although seven morphometric characteristics were measured, only four (beak length, culmen length, beak width, and tarsal length) were used for the analyses, as it was not possible to measure the other three (culmen depth, wing length, and tail length) in many samples because of damage to the beaks, wings, or tails in the specimens. All measurements, except for the culmen length, were taken by C-M Hung. Culmen length was measured by S-H Li. Student’s t-tests were used to determine whether the morphometric characteristics differed between the two subspecies; male and female samples were compared separately.

Results

Genetic diversity of B. thoracicus and B. sonorivox

We sequenced 21,150 base pairs of 31 exon loci in total. The level of nucleotide diversity (π and Θ_w) of B. thoracicus was higher than that of B. sonorivox (Additional file 1: Table S3). The values of these summary statistics all suggest that the effective population size of B. thoracicus should be larger than that of B. sonorivox. The genetic differentiation (F_ST) between the Bamboo Partridges was high based on the 31 nuclear loci examined (mean = 0.432; range: 0.025–0.889). Results of the multilocus HKA test showed that the polymorphisms observed in all loci did not significantly deviate from the neutral expectation (p > 0.05). The estimated substitution rates of the nuclear loci ranged from 1.94 × 10 ^‒9 to 9.31 × 10⁻⁹ substitutions/site/year (mean = 5.19 × 10⁻⁹ substitutions/site/year). The numbers of substitutions per locus per year was shown in Additional file 1: Table S3 with other summary statistics for all 31 nuclear loci.

Historical demography since the split of the two diverging Bamboo Partridges

The ABC model suggested that the early gene flow model was favoured (posterior probability = 0.53; Fig. 2) by the 31 nuclear markers, supporting the idea that post-divergence gene flow only occurred within the first two tenths of the divergence history. Our results cannot fully reject the late gene flow model (posterior probability = 0.30); the Bayes factor of the early interruption vs. late interruption is only 1.78. Therefore, the late interruption model also received some support from our polymorphism data.

Both early and late gene flow models generated similar effective population sizes for the two Bamboo Partridges and their common ancestor (median N_Bt = 5.12‒5.89 × 10⁵, median N_Bs = 2.47‒2.67 × 10⁵, and median N_a = 9.63‒9.76 × 10⁵; Table 1). The time of splitting between the two Bamboo Partridges was estimated to be approximately 1.21 (early gene flow model) to 1.43 (late gene flow model) million years ago. For the most likely early gene flow model, the post-divergence gene flow between the two diverging Bamboo Partridge lineages ceased approximately 1.12 × 10⁶ years ago (Table 1), while the cessation was approximately 0.63 × 10⁶ years ago for the late gene flow model. In addition, a low migration rate was suggested by both models: for the early gene flow model, the medians of M_Bs and M_Bt were 0.0049 and 0.0051, respectively; for the late gene flow model, they were 0.0013 and 0.0043, respectively (Table 1).

Table 1 Population parameters of the two models estimated by msABC

Ecological niches of the two Bamboo Partridges

The ENMs produced by Maxent showed that the predicted distributions of both Bamboo Partridges closely matched their actual occurrences (Additional file 2: Figures S4, S5). The mean value of the AUC was similar for the two species (0.805 for B. thoracicus and 0.783 for B. sonorivox), suggesting that the ENM results were reliable.

According to the percent contribution, the ecological niches of the two Bamboo Partridges seem to be determined by different climatic variables: the range of B. thoracicus is mainly determined by Bio 2 (36.5%) and Bio 15 (24.1%), while the range of B. sonorivox is mainly influenced by Bio 9 (22.0%), Bio 10 (20.2%), Bio 12 (21.8%), and Bio 7 (17.4%). However, the permutation importance results were not consistent with those of the percent contribution, the range of B. thoracicus is mainly influenced by Bio 9 (24.9%) and Bio 7 (15.2%), and the range of B. sonorivox is more influenced by Bio 10 (28.7%) and Bio 7 (18.4%).

Ecological niche differentiation between the two Bamboo Partridges

The actual Schoener’s D (0.047) and similarity index I (0.209) suggested that the ecological niches of B. thoracicus and B. sonorivox were highly differentiated. In contrast, the background test suggested that ecological niches of the two Bamboo Partridges were more conservative than expected from the stochastically simulated ecological niche in the geographic range of B. thoracicus (99% confidence intervals of the simulated values: Schoener’s D = 0.005–0.035, I = 0.033–0.157, p < 0.01) and B. sonorivox (99% confidence intervals of the simulated values: D = 0.000–0.011, I = 0.001–0.760, p < 0.01) (Fig. 3). The results of the background test implied that the ecological niches of the two Bamboo Partridges were not significantly different, and the set of important variables that differed between the two ENMs were likely caused by high climatic habitat heterogeneity within the range of two Bamboo Partridges.

Fig. 3

The background similarity test shows that B. thoracicus and B. sonorivox have significantly conserved ecological niches. In both a and b, the y-axis means the frequency of the similarity index (D and I) from 200 random simulations. D and I reflect the extent of niche overlap. The red histogram shows the distribution of similarity index, comparing the actual ecological niche of B. sonorivox with the simulated ecological niche of B. thoracicus. The grey histogram shows the distribution of similarity index in a reversal scenario.

Potential paleo-distribution of B. thoracicus and B. sonorivox during the late Pleistocene

By extrapolating ENMs of B. thoracicus and B. sonorivox to the paleoclimatic data of East Asia during the LGM, we found that, compared to its current range, the potential distribution of B. thoracicus during LGM was reduced and its potential habitable region on Taiwan Island was smaller and was disconnected from the other major suitable regions in southwestern China and Northern Indochina (Fig. 4; Additional file 2: Figure S4). Our results indicate that B. sonorivox could expand its range in the LGM (distributed on both the mainland and Taiwan Island) compared to its current range (Fig. 4). Because 75% of the suitable region of B. sonorivox overlapped with the potential range of B. thoracicus in Taiwan (Fig. 4), the ancestors of the two Bamboo Partridges were likely to have had secondary contact during the last glacial period.

Fig. 4

The potential distribution area of B. thoracicus and B. sonorivox in the LGM. Dark green area shows the potential secondary contact region of the two species in the LGM. Brown and yellow areas show the potential suitable distribution zones of B. thoracicus and B. sonorivox in the LGM, respectively. The grey color region shows the unsuitable area for both species in the study region.

Morphological differentiation between B. thoracicus and B. sonorivox

The only significant difference in morphological characteristics we observed was in beak length (p ≤ 0.01; Table 2). The average beak length of B. thoracicus (♂) was 20.74 mm, which is smaller than that of B. sonorivox (average: 22.02 mm, ♂). For the female specimens, the average beak length of B. thoracicus was also smaller than that of B. sonorivox (19.77 vs. 21.76 mm). For the other three morphological characteristics, beak width, culmen length, and tarsus length, there was no significant difference between B. thoracicus and B. sonorivox of the same gender.

Table 2 Summary of the morphological differences between B. thoracicus and B. sonorivox

Discussion

In the current project, we employed 31 exon loci with the ABC model to estimate the duration of the post-divergence gene flow, and used ENM to predict the potential distribution of the two Bamboo Partridges in the LGM. The results suggested that gene flow occurred between the two Bamboo Partridge species during the divergence process, and that they might have shared the same distribution during the LGM.

Speciation with gene flow

Our genetic data suggested that there was a low-level gene flow (M for both directions below 0.01; Table 1) between ancestral lineages of the Chinese and Taiwan Bamboo Partridges in both the early and late stages of their speciation (Fig. 2; Table 1). Even though the magnitude of gene flow was low, the evidence did not support the isolation model (Fig. 2). The Chinese and Taiwan Bamboo Partridge populations had gene flow with each other during their divergence process (Fig. 2; Table 1), which is consistent with the case of Taiwan Hwamei (Leucodioptron taewanus) and Hwamei (L. canorum canorum) (Li et al. 2010). The Taiwan Hwamei, distributed on Taiwan Island, also had gene flow with its mainland sister species, Hwamei, during their divergence process (Li et al. 2010). The two Bamboo Partridge species began to diverge from each other during the early Pleistocene, and the gene flow ceased in the mid-Pleistocene (Table 1). The fluctuations in the sea level in the Taiwan Strait during early late Pleistocene may have led to changes in their distribution alternating between sympatry and allopatry (Shi et al. 2006), which may promote their divergence (He et al. 2019). The complex fluctuations in sea level could lead to the cycles of mixing-isolation-mixing, which further contribute to the speciation (He et al. 2019). As shown in our ENM analysis, both species shared the same distribution during the glacial period (Fig. 4), followed by an allopatric distribution during the inter-glacial period. Multiple glaciation cycles during the last million years may have led to the divergence of the two Bamboo Partridge species, during the cycles of allopatric and sympatric distribution resulting from fluctuations in sea level.

The low level of post-divergence gene flow could be explained by two alternative hypotheses. At first, it could have been initiated by a mutation affecting plumage colouration that resulted in effective self-matching assortative mating. The rise of self-matching plumage traits could have effectively reduced the level of gene flow between ancestral Bamboo Partridges, even in sympatry. However, a caveat of such a scenario is that the Bamboo Partridges with different plumage colouration would compete for the same ecological niche. Observations carried out in captivity and field survey data suggest that the two species utilize similar food types (Kobayashi 1976; Zheng 2015). Unless their different beak lengths allowed partitioning of their realized niches effectively in sympatry, their long-term coexistence would be highly unlikely. Alternatively, it is possible that the two Bamboo Partridges diverged from a highly structured ancestral species. It has been shown that the genetic signature of speciation with gene flow between two populations is indistinguishable from allopatric speciation from a subdivided ancestral population (Yang et al. 2017). If two species diverged from partially geographically differentiated populations or a structured population, the ancestral population should, in theory, have a much larger effective size than that of the extant population for its deep structure (Yang et al. 2017). For the Chinese and Taiwan Bamboo Partridges, the effective population size of the ancestral species is estimated to be much larger than that of both extant species (Table 1). Therefore, it is possible that the two Bamboo Partridges diversified from a geographically structured ancestral species. The deep structure of ancestral species greatly diminished the level of gene flow between ancestral populations of the two species at this stage. Mating characteristics and other trophic-related functional traits diverged later on.

Ecological niche conservatism in the two Bamboo Partridges

Accumulating recent data suggest that divergence of ecological niches could be the force driving species splitting or maintenance of species boundaries between incipient species (Shaner et al. 2015; Wang et al. 2019; Germain et al. 2020). However, fundamental niches, such as those pertaining to the climatic, are usually considered to be highly conserved (Peterson et al. 1999), especially for recently diverged lineages (Peterson 2011). It has been suggested that niche conservation could promote speciation when ancestral lineages are isolated by unsuitable habitats (Wiens 2004; Gutiérrez-Ortega et al. 2020). In other words, niche conservatism could be a critical component in driving allopatric speciation.

However, the presence of gene flow during the early phase of divergence does not fit with the classical allopatric speciation model in the two Bamboo Partridges (Mayr 1942, 1963). The potential paleo-distribution of the two Bamboo Partridges in the LGM (Fig. 4) also suggested a limited role of climatic, or fundamental niche divergence in the speciation of the B. thoracicus and B. sonorivox. Our projections of the ENMs to the LGM suggested that the ranges of the ancestors of the two species could have overlapped (Fig. 4). The conserved ecological niche of the two Bamboo Partridges may provide numerous potential opportunities for secondary contact during the dry and cool glacial periods in the Pleistocene. The conserved ecological niche may contribute to the allopatric speciation based on theory; however, it may also contribute to the speciation with gene flow during secondary contact.

Potential evolutionary processes that contribute to reproductive isolation of the two Bamboo Partridges

The Chinese and Taiwan Bamboo Partridges have reproductive isolation based on the observation in the zoo. In Japan, no hybrid has been reported, even where the two species are feed in the same cage, suggesting complete reproductive isolation (Kobayashi 1976). Both environmental selection and sexual selection, can promote species divergence in the face of gene flow (van Doorn et al. 2009; Richards et al. 2019; Stankowski et al. 2019). Environmental selection can reduce the incidence of mating between individuals from different populations or reduce the fitness of hybrids through ecological niche variation between parental lineages. However, the conserved ecological niche of the two Bamboo Partridges (Fig. 3) suggested the divergent environmental selection may have less power on driving speciation, because the environmental traits of their habitat may be similar.

Assortative mating can rapidly drive morphological differences at the initial stage of speciation (Higashi et al. 1999; Boughman 2001; Castillo et al. 2015; Lindsay et al. 2019). However, some theoretical and empirical evidence suggest that sexual selection can only work with environmental selection or genetic drift to achieve complete speciation (van Doorn et al. 2009; Lindsay et al. 2019). Without the participation of environmental selection, it would be difficult for assortative mating to be maintained in the same ecological niches: the ancestral species could hybridize with each other within an ecological niche unless the intensity of sexual selection was extremely strong (Arnegard and Kondrashov 2004; van Doorn et al. 2009; McLean et al. 2020). Therefore, sexual selection must interact with disruptive ecological selection in the speciation process (van Doorn et al. 2009; McLean et al. 2020). The beak is a critical functional trait in birds for capturing and processing food (Cooney et al. 2017) and is therefore highly associated with their trophic or realized niche (Pigot et al. 2020). Our finding that the beak lengths of B. thoracicus and B. sonorivox are significantly different (Table 2) implies that the two species could utilize different realised niches associated with their food resources. Therefore, assortative mating might work synergistically with environmental selection to facilitate speciation of the two Bamboo Partridges. The speciation gene uncovered by the genetic data may help us to further explore the potential drivers of speciation between the two Bamboo Partridges (Yang et al. 2017; Wang et al. 2020). If the speciation gene features linked with the mating behaviour or male phenotype, assortative mating may drive speciation (Ting et al. 2001; Wang et al. 2020).

Conclusions

The current project provided evidence of a conserved ecological niche, different beak lengths, and post-divergence gene flow, which suggested that gene flow occurred between allopatric Chinese and Taiwan Bamboo Partridges during the speciation process. The evidences supported speciation with gene flow model for the allopatric species. Sexual selection may work synergistically with environmental selection to drive reproductive isolation. The results of the current project will improve our understanding of the avian speciation process.