Incipient speciation with gene flow on a continental island: Species delimitation of the Hainan Hwamei (Leucodioptron canorum owstoni, Passeriformes, Aves)

Accepted Manuscript Incipient speciation with gene flow on a continental island: Species delimitation of the Hainan Hwamei (Leucodioptron canorum owstoni, Passeriformes, Aves) Ning Wang, Bin Liang, Jichao Wang, Chia-Fen Yeh, Yang Liu, Yanlin Liu, Wei Liang, Cheng-Te Yao, Shou-Hsien Li PII: DOI: Reference:

S1055-7903(16)30111-7 http://dx.doi.org/10.1016/j.ympev.2016.05.022 YMPEV 5531

To appear in:

Molecular Phylogenetics and Evolution

Received Date: Revised Date: Accepted Date:

7 November 2015 26 April 2016 20 May 2016

Please cite this article as: Wang, N., Liang, B., Wang, J., Yeh, C-F., Liu, Y., Liu, Y., Liang, W., Yao, C-T., Li, SH., Incipient speciation with gene flow on a continental island: Species delimitation of the Hainan Hwamei (Leucodioptron canorum owstoni, Passeriformes, Aves), Molecular Phylogenetics and Evolution (2016), doi: http:// dx.doi.org/10.1016/j.ympev.2016.05.022

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Incipient speciation with gene flow on a continental island: Species delimitation of the Hainan Hwamei (Leucodioptron canorum owstoni, Passeriformes, Aves) Ning Wanga1*, Bin Lianga1, Jichao Wanga, Chia-Fen Yehb, Yang Liuc, Yanlin Liud, Wei Lianga, Cheng-Te Yaoe, Shou-Hsien Lib*

1

a

Ministry of Education Key Laboratory for Tropical Plant and Animal Ecology,

College of Life Sciences, Hainan Normal University, Haikou 571158, China. b

Department of Life Science, National Taiwan Normal University, Taipei 11677,

Taiwan. c

State Key Laboratory of BioControl, College of Ecology and Evolution/School of

Life Sciences, Sun Yat-sen University, Guangzhou 510275, China. d

Institute of Forestry Ecology, Environment and Protection, Chinese Academy of

Forestry, Beijing 100091, China e

High-altitude Experimental Station, Taiwan Endemic Species Research Institute,

Chi-chi 55244, Taiwan.

Corresponding authors: Ning Wang: 99th south Longkun Rd., Haikou 571158, Hainan, China. E-mail address: [email protected]. Shou-Hsien Li: 88 Sec. 4, Tingchou Rd., Taipei 11677, Taiwan. E-mail address: [email protected].

1

Ning Wang and Bin Liang contribute equally to this work.

2

ABSTRACT Because of their isolation, continental islands (e.g., Madagascar) are often thought of as ideal systems to study allopatric speciation. However, many such islands have been connected intermittently to their neighboring continent during recent periods of glaciation, which may cause frequent contact between the diverging populations on the island and continent. As a result, the speciation processes on continental islands may not meet the prerequisites for strictly allopatric speciation. We used multiple lines of evidence to reevaluate the taxonomic status of the Hainan Hwamei (Leucodioptron canorum owstoni), which is endemic to Hainan, the largest continental island in the South China Sea. Our analysis of mitochondrial DNA and twelve nuclear loci suggests that the Hainan Hwamei can be regarded as an independent species (L. owstoni); the morphological traits of the Hainan Hwamei also showed significant divergence from those of their mainland sister taxon, the Chinese Hwamei (L. canorum). We also inferred the divergence history of the Hainan and Chinese Hwamei to see whether their divergence was consistent with a strictly allopatric model. Our results suggest that the two Hwameis split only 0.2 million years ago with limited asymmetrical post-divergence gene flow. This implies that the Hainan Hwamei is an incipient species and that speciation occurred through ecologically divergent selection and/or assortative mating rather than a strictly allopatric process.

Keywords: species delimitation, incipient speciation, post-divergence gene flow, continental island.

3

1. Introduction Speciation, the process creating biodiversity, has been one of the central issues of biology (Coyne and Orr, 2004) since Darwin published his “On the Origin of Species” (Darwin, 1859). Allopatric isolation, the complete cessation of gene flow between populations isolated by geological events, as the common process by which new species arise (Lande, 1980, Mayr, 1942), could allow reproductive isolation to arise as a byproduct of genetic divergence through local adaptation and/or genetic drift (Hoskin et al., 2005; Mayr, 1942). However, the efficiency of a geological barrier in isolating a population also correlates with the stability of its existence and the mobility of the organism (Campbell and Reece, 2002; Gaston, 2003). Although continental islands (such as the Aegean Archipelago, Madagascar) are usually treated as systems to study the effects of geographical isolation on speciation (e.g., Comes et al., 2008; Whittaker and Fernández-Palacios, 2007), the geographical barrier between these islands and the nearby continent may be inconsistently present, as water bodies dry out periodically due to the depression of the sea level during glacial periods (e.g., Lambeck and Chappell, 2001) allowing intermittent genetic exchange among the once-isolated populations. As a consequence, speciation with the postdivergence gene flow (such as parapatric speciation, e.g., Barraclough and Vogler, 2000) becomes a viable alternative model for the formation of continental island endemics. The homogenizing effect of gene flow has been thought to be too strong to overcome the forces of divergence (e.g., Slatkin, 1987; but see Garant et al., 2005, Gavrilets and Vose, 2005), but post-divergence gene flow now seems more common in speciation processes than previously assumed (e.g., Bank et al., 2012; Li et al., 2010a; Niemiller et al., 2008;

4

Won and Hey, 2005; Yeung et al., 2011). To test these alternative hypotheses of continental island speciation, we used evidence from multiple sources to evaluate the taxonomic status and the level of postdivergence gene flow of an endemic resident passerine in Hainan, the Hainan Hwamei (Leucodioptron canorum owstoni), which occupies the margins of lowland secondary evergreen woodlands (MacKinnon and Phillipps, 2000). As the largest continental island in the South China Sea, Hainan possesses a large amount of subtropical and tropical forest and hosts a variety of endemic biota (Xing et al. 1995, Zheng, 2011). Hainan is separated from the Leizhou Peninsula of the Asian mainland by the shallow and narrow Qiongzhou Strait with an average width of 30 km and a maximum depth of ~120m (Shi et al., 2002). During the Pleistocene, the global sea level fluctuated with each glacial cycle (Williams et al., 1998), causing a connection-disconnection oscillation of this mainland-island system (Yan, 2006). Consequently, Hainan has been connected with the nearby continent for approximately 50% of the last two hundred and fifty thousand years (Voris, 2000). Therefore, isolation of the gene pools on the continent and Hainan should have been restricted largely to interglacial periods with the potential contact when a land bridge formed in glacial periods (Lambeck and Chappell, 2001; Voris, 2000) . Since first described by Rothschild in 1903, the Hainan Hwamei (L. c. owstoni) has conventionally been considered as a subspecies of the Chinese Hwamei (L. canorum) (e.g., Clements et al., 2015; MacKinnon and Phillipps, 2000). The nominate subspecies of Chinese Hwamei, L. c. canorum, is widely distributed in continental East Asia, ranging from northern Indochina to central and southern China (Fig. 1). Leucodioptron c. owstoni is described as having considerably paler plumage than L. c. canorum: its upper-side is

5

distinctively more olive and less yellowish; its ear cover is duller and darker (Rothschild, 1903). A third taxon in genus Leucodipotron is the Taiwan Hwamei (L. taewanum, Swinhoe, 1859) present in the lowlands of Taiwan. Leucodioptron taewanum used also to be treated as a subspecies of L. canorum (e.g., Howard and Moore, 1994), but has been granted full species status based on both morphological and molecular evidence (Li et al., 2006). Although an analysis based on mtDNA sequences suggested that L. c. owstoni forms a monophyletic clade that separated from L. c. canorum around 600 thousand years ago (Kya), the limited sample size of the previous study (Li et al., 2006) rendered a conclusion about its taxonomic status uncertain. It has been suggested that falling sea levels during repeated glacial periods could have facilitated gene exchange between the diverging L. canorum and L. taewanum (Li et al., 2010a), which is consistent with the incomplete reproductive isolation of the two species (Li et al., 2010b; Tu and Severinghaus, 2004). Given the even shorter distance between Hainan and the adjacent continent, and the closer phylogenetic relationship between L. c. canorum and L. c. owstoni (Li et al., 2006), it is possible that reproductive isolation between them is also incomplete. Therefore L. c. canorum and L. c. owstoni should serve an ideal system to examine whether post-divergence gene flow occurred between a recently diverged pair of continent and continental-island taxa. Here, we used twelve nuclear regions, together with the mitochondrial cytochrome b (CYTB) gene to evaluate the taxonomic status of the Hwamei complex based on both the gene tree and species tree (Heled and Drummond, 2010), population structure and coalescent-based species delimitation (Yang and Rannala, 2010) methods. Moreover, morphometric characters measured from skin specimens were compared to detect

6

morphological divergence between the Hainan Hwamei and their mainland relatives. We also applied the isolation-with-migration model (Hey and Nielsen, 2007) to estimate their divergence time and the level of post-divergence gene flow between them (allowing for the presence of post-divergence gene flow to avoid underestimating divergence time, Leaché et al., 2014). Consistent genetic and morphological analyses suggested that L. c. owstoni could be treated as a full species that recently diverged. This study not only provides a case of a non-allopatric speciation, but also illustrates the prevalence of postdivergence gene flow in shaping the species-divergence processes in continental-island systems.

2. Material and methods

2.1. Sample preparation Blood, feather pulp (i.e., the calamus section of the feather) or muscle samples from 13 L. c. owstoni individuals were collected in Diaoluo Moutain (N=10), Wanling (N=2), and Jianfengling (N=1) in Hainan Island and used to extract the gross genomic DNA. Blood, liver or muscle samples of 79 L. c. canorum individuals had previously been collected (e.g., Li et al. 2010a) from part of its entire range, while 30 L. taewanum individuals were caught in central Taiwan (Fig. 1, Appendix A, Table A.1). We used one sample of Garrulax monileger schmacheri (Gms1, from Wanling, Hainan) as the outgroup when conducting phylogenetic analyses. DNA was extracted using the whole genome DNA extraction kit (Vicband Life Sciences company [HK] Limited), suspended in Elution Buffer, and stored at −20 °C.

7

In total, 12 nuclear markers were designed from the homologous nuclear regions between chicken (Gallus gallus, Hillier et al., 2004) and zebra finch (Taeniopygia guttata, Warren et al., 2010) (Table A.2). The entire mitochondrial CYTB gene was also amplified via polymerase chain reactions (PCRs) by using the L13653/H14296 and L14192/H14853 primers (Li et al., 2006). We used 40 l PCR reaction volumes, including ~200 ng DNA, 10×PCR Buffer, 0.2 mM of each dNTP, 0.15 M of each primer, 0.75 U Taq DNA polymerase (Amersham Biosciences) and 1.5 mM MgCl2. The PCR profile for CYTB was denaturation at 94°C for 5 min, 35 cycles of 94°C for 30 s, 55°C for 30 s and 72°C for 1 min, followed by final extension at 72°C for 10 min. All nuclear loci were amplified using a single touchdown PCR protocol (The annealing temperature was from 60°C to 50°C, being decreased by 0.5°C per cycle followed by 20 additional cycles at 50°C). All PCR products were examined for size on 1% agarose gel and purified using PCR purification Kit (QIAGEN). Samples of L. c. owstoni were sequenced in both directions using ABI BigDye® Terminator v.3.1 on an ABI PrismTM 3100-Avant genetic analyzer in Shenzhen Genomic Institute. Samples of L. c. canorum and L. taewanum were sequenced for short reads and mapped to the reference sequences from a single individual of L. taewanum generated by Sanger sequencing (LifeScopeTM Genomic Analysis Software from Life Technologies), mainly following the protocol of Shaner et al. (2015). The novel sequences collected in this study have been deposited in Genbank (Table A.1).

2.2. Population genetics analyses The nuclear DNA genotypes were reconstructed using the PHASE (Stephens et al.,

8

2001) algorithm under a maximum likelihood method implemented in DnaSP v5 (Librado and Rozas, 2009). Alleles with probabilities lower than 60% (Chu et al., 2013; Harrigan et al., 2008; Hung et al., 2014) were excluded from the subsequent analyses when the phased data were required (e.g., Structure and IMa analyses). The number of segregating sites, the number of alleles, nucleotide diversity (π) and Watterson’s θ (θw) were calculated. The nucleotide distances (Dxy; Nei, 1987) between L. c. canorum and L. c. owstoni and between L. c. canorum and L. taewanum were also calculated. We tested the three Hwamei for selective neutrality using the Hudson–Kreitman–Aguade (HKA) test (Hudson et al., 1987) and Tajima’s D statistics (Tajima, 1989). All the above analyses were conducted in DnaSP v5. We also looked for evidence of recombination using the Phi test (Bruen et al., 2006) implemented in the program SplitsTree 4.13.1 (Huson and Bryant, 2006), which is designed to distinguish recombination from recurrent mutation.

2.3. Phylogenetic analyses- gene tree approach After aligning the sequences with Mafft 6.717 (Katoh et al., 2009), we made a concatenated dataset using combine-0.9 (written by Edward L Braun) by first including then excluding the CYTB gene with all the un-phased nuclear sequences (i.e., 13 loci and 12 loci). In order to reduce the missing data effect, 39 CYTB gene sequences that were collected previously (Table A.3) were not used in the concatenated dataset, though they were included in the subsequent analyses of the single CYTB gene. Using the two sets of concatenated matrices, we conducted the partitioned maximum likelihood (ML, assigned each locus with independent substitution model) analyses in RAxML 7.2.6 (Stamatakis,

9

2006) with the GTR+ model and ten randomized starting trees. To test the robustness of phylogenetic inference, partitioned ML bootstrap analyses with 500 replicates using the same model were conducted in RAxML. As for the single CYTB gene, we reconstructed the phylogenetic relationships among individuals by using Bayesian inference (BI) in MrBayes 3.1.2 (Ronquist and Huelsenbeck, 2003) under the best model determined by the Akaike information criterion (AIC) in Modeltest 3.7 (Posada and Crandall, 1998). Two independent Markov chain Monte Carlo (MCMC) runs were performed, each containing one cold chain and three heated chains, for 100 million generations and sampled every 10,000 generations. The convergence between runs was examined using the average standard deviation of split frequencies (i.e., <0.01) and the initial 25% of each run was discarded as burn-in. Using the CYTB gene, individual-based ML trees were also built in PAUP*4.0b10 (Swofford, 2002) (with GTR++I model and ten replications) and RAxML (using GTR+ model and ten random starting trees). To incorporate allelic variation from multiple nuclear markers (Joly and Bruneau, 2006), an individual-based, multi-locus genetic network (e.g., Barrett and Freudenstein, 2011; Dong et al., 2014) was constructed by using all 12 nuclear loci after excluding individuals for which more than four loci were missing. The uncorrected p distances for all the loci were calculated first using PAUP*, and were then combined into one distance matrix using POFAD 1.07 (Joly and Bruneau, 2006). With this, a genetic network of Hwamei individuals was constructed using a NeighborNet algorithm (Bryant and Moulton, 2004) in SplitsTree.

10

2.4. Phylogenetic analyses- species tree approach We treated all three taxa in genus Leucodioptron as operational taxonomic units and jointly estimated multiple gene trees embedded in a shared species tree using a Bayesian MCMC method implemented in the program *BEAST 1.8.1 (Drummond and Rambaut, 2007). Analyses were conducted for both the 12 nuclear genes and all 13 markers. We assigned samples collected from the same taxon to the same population. The optimal or the closest model of nucleotide substitution was assigned to each locus based on the AIC in Modeltest. Uncorrelated lognormal relaxed clock models were used for each locus (Drummond et al., 2006). When the mitochondrial gene was included in the species tree analysis, we set a clock rate prior of 1.035% (0.01035/site/million years) for CYTB, following a previous suggestion for passerines (Weir and Schluter, 2008). When only nuclear genes were used in the analyses, we set a clock rate prior of 0.32% (0.0032/site/million years) for locus ADIPOR1; this rate was calculated from the CYTB rate using the same protocol as in Li et al. (2010a). The prior distribution of clock rate followed a lognormal distribution. We used a Yule process as the species tree prior, and a piecewise linear and constant root as the population size model. Two independent runs with 500 million generations were conducted (sampling every 50,000 generations and discarding the first 25% as burn-in) in *BEAST. The convergence of these two chains was assessed in the program Tracer 1.4.1 (Rambaut and Drummond, 2007).

2.5. Population structure analyses Pairwise differentiation between all three Hwamei subspecies/species was estimated, based on both CYTB and the un-phased concatenated nuclear genes, by

11

computing pairwise FST (Hudson et al., 1992) using the pairwise number of nucleotide differences with the Weir and Cockerham (1984) estimator in Arlequin 3.11 (Excoffier et al., 2005). The significance of the associated statistics was assessed using 10,000 nonparametric permutation procedures (Excoffier et al., 1992). We tested the genetic structure of the three Hwamei taxa using a Bayesian clustering method implemented in Structure 2.3.4 (Falush et al., 2003; Pritchard et al., 2000) based on the 12 phased nuclear genes. After data format conversion in PGDSpider (Lischer and Excoffier, 2012), we conducted the Structure analysis by using an admixture model with correlated allele frequencies among populations, and performing 3×10 5 MCMC steps after a 3×105-step burn-in period. In order to detect potential cryptic genetic structure, ten independent runs for each K-value (K=1-7) were conducted for the entire dataset. We used Structure Harvester online program (Earl and vonHoldt, 2012) to identify the most likely number of genetic clusters based on the statistics (Evanno et al., 2005). The final results of the bar plot for individual memberships were drawn with a cluster visualization program Distruct (Rosenberg, 2004).

2.6. Coalescent-based species delimitation Joint Bayesian species delimitation and species tree estimation was conducted using BP&P (Yang, 2015) under the multispecies coalescent model in a Bayesian framework. It can account for incomplete lineage sorting due to ancestral polymorphism and gene treespecies tree conflicts (Yang and Rannala, 2010; Yang and Rannala, 2014). We used the un-phased nuclear sequences in the BP&P analyses. Both diffuse (2, 200) and informative (10, 1000) gamma priors, with mean 2/200 (10/1000) = 0.01 (one difference

12

per 100bp), were used for θ (θ= 4Neµ, Ne: effective population size, µ: substitution rate per site per generation). We also assigned a diffuse gamma prior (2, 600) to the ancestral root age τ0, while the other τ were assigned with the Dirichlet prior (Yang and Rannala, 2010). The prior mean for τ0 (τ0 = tµ, t is the root age in generations) was approximated to the divergence time estimated from the *BEAST analysis (using 12 nuclear loci) with an assumed generation time of 2.5 years (Li et al., 2010a). For µ, we used 2.3×109

/site/year based on the average substitution rate of nuclear loci relative to that of CYTB.

We conducted the BP&P analyses for 550,000 generations, sampling every five generations, and discarded the first 50,000 generations as burn-in. Each prior set was run twice to confirm consistency between runs.

2.7. Divergent demography between L. c. canorum and L. c. owstoni The effective population size parameters (θ= 4Neµ, µ: substitution rate per gene per generation), population divergence times (t= Tµ, T: time since divergence in generations) and directional migration rates (m= m/µ, m: the proportion of immigrant gene copies in a population per generation) between L. c. canorum and L. c. owstoni were estimated based on the phased nuclear data with/without CYTB using the coalescent-based MCMC method implemented in IMa (Hey and Nielsen, 2007). To convert the scaled parameters, we used a substitution rate of 1.035×10-8/site/year for CYTB with a confidence interval (CI) between 0.6-1.505×10-8/site/year as in previous studies (Weir and Schluter, 2008). We also assigned mutation rates to nuclear loci whenever they could be estimated by comparison with the CYTB substitution rate based on the average genetic distances between each subspecies and the outgroup in the current study (Table 2). We assumed the

13

generation time as 2.5 years for Hwamei (Li et al., 2010a). After deleting the first one million steps as burn-in, two independent MCMC runs of at least 20 million generations were conducted under a geometric heating scheme (-fg –n40 –g1 0.8 –g2 0.9). The effective size values (ESS > 300) and the plots of trend lines were examined to confirm the convergence of parameter estimates.

2.8. Morphometric analyses Bin Liang recorded nine morphometric characteristics (Table A.4) of 58 male Hwamei specimens collected from southern China, including Guangdong (N=14), Guangxi (N=20), Yunnan (N=5), and Hainan (N=19) provinces, and archived in the Museum of Biology, Sun Yat-sen University, and the National Zoological Museum of China, Institute of Zoology, Chinese Academy of Sciences (NZMC, IOZ, CAS). Seven characteristics [beak length (BeL), tarsal length (TaL), wing length (WL), tail length (TL), averaged length of both eyebrows (aEL, from the beginning of eye-ring to the end of eye line), original body length (oBL) and weight] were compared between L. c. owstoni and L. c. canorum from each province and as a whole, and Student’s t calculated. We also applied linear discriminant analysis, which can determine how well distinct two or more groups are based on several characteristics of those groups (Corstanje et al., 2009, Liao et al., 2014), to all variables except for Weight and oBL (which contained too many missing data). All statistical analyses were conducted in JMP 12 (SAS Institute Inc., Cary, NC, USA).

3. Results

14

In total, we obtained 1,143 bp of the CYTB gene and 6,475 bp of the nuclear sequences. No internal stop codons or indels were found in CYTB, suggesting the absence of a nuclear copy of the mitochondrial gene (e.g., Sorenson and Fleischer, 1996), in the sequencing results. We found more alleles (haplotypes) in the CYTB gene than in any nuclear sequence (Table 1). For most of the nuclear loci, L. c. owstoni exhibited higher nucleotide diversity and genetic diversity (θw) than L. taewanum (Table 1), although L. c. canorum always showed the highest diversity (This could be due to sampling bias, as more individuals were sampled from this taxon). Moreover, the nucleotide differences between L. c. owstoni and L. c. canorum were higher in five out of the 13 loci than those between L. c. canorum and L. taewanum, although the opposite pattern was shown in seven out of the other eight loci (Table 1). Based on Tajima’s D statistics, the null hypothesis of evolutionary neutrality could only be rejected at a significance level of 0.05 for SDC2, and not for any of the other genes. When using the HKA test, no locus exhibited evidence of selection. No evidence of recombination within a locus was found in the Phi test (Table 1).

3.1. Incongruence among the topology of phylogenetic trees Partitioned ML analyses based on the concatenated dataset (including both nuclear and mtDNA) suggested that L. taewanum forms a monophyletic clade (bootstrap support 76%, Fig. 2A) basal to the other Hwameis. However, L. c. canorum formed a paraphyletic group with respect to the monophyletic L. c. owstoni, which only obtained marginal support (18%) in the RAxML analyses (Fig. 2A). When excluding CYTB, the

15

partitioned ML analyses based on the concatenated 12 nuclear DNA suggested that neither L. c. owstoni or L. c. canorum was monophyletic, although L. taewanum still existed as a monophyletic clade with relatively high support (71%, Fig. 2B). However, L. taewanum nested within L. c. canorum rather than forming the basal lineage [Fig. 2B, although its basal position was supported in the bootstrap consensus tree (Appendix B)]. In order to detect the potential impact of individual nuclear gene on such rearrangement, we conducted gene-jackknifing analyses (Hackett et al., 2008) by excluding one gene at a time and running ML tree in RAxML with the rest of the nuclear dataset. We found that four loci (F1N9H4, ABCC5, MOCS1 and HPGDS) might contribute to this topological rearrangement. When excluding single genes or all four of them, L. taewanum again formed the basal clade (Appendix B), suggesting that its basal position is not driven exclusively by CYTB. ML analysis in both RAxML and PAUP* (not shown) based on the single CYTB gene presented a clear three-clade phylogeny (Fig. 3A), in which L. taewanum split first, with a later separation (bootstrap support 41%) between L. c. canorum and L. c. owstoni. The same pattern was found in the BI analysis. Moreover, the inferred multi-locus nuclear network supported the genetic distinctness of L. taewanum, and exhibited a strong structure of L. c. owstoni separating from the mainland group (Fig. 3B). The *BEAST species tree based on either the 13 loci or the 12 nuclear loci strongly supported (posterior probability=1) a sister relationship between exclusive monophyletic L. c. canorum and L. c. owstoni, which further group with L. taewanum (Fig. 3C).

3.2. Genetic differentiation and species delimitation

16

Genetically, all three Hwamei taxa are highly divergent from each other. Using CYTB data, the pairwise FST between L. c. canorum and L. c. owstoni was significantly larger than zero (FST= 0.813, P= 0.000), which was comparable with those between L. c. canorum and L. taewanum (FST= 0.901, P= 0.000) and between L. c. owstoni and L. taewanum (FST= 0.927, P= 0.000). A similar pattern was also found when using concatenated nuclear DNA, with the pairwise FST between L. c. canorum and L. c. owstoni (FST= 0.863) even higher than that between L. c. canorum and L. taewanum (FST= 0.835, Table 3). The genetic distinctiveness of the three Hwamei taxa is also supported by the results of the Bayesian clustering analysis. The results of the Structure analysis based on the nuclear dataset revealed three genetic clusters (Fig. 3D), each representing a specific Hwamei taxon in all 10 replicates (The highest peak of K occurred at K= 3, Fig. A.1). Our results also suggest that no Hainan samples used in the current study should be regarded as a hybrid or admixed individuals (assignment probability > 0.97 for all individuals). It has been suggested that high posterior probabilities (> 0.95) in the BP&P analyses indicate greater statistical support for resolving a certain node into multiple species (Leaché and Fujita, 2010). In the current study, the posterior probability for delimitating L. c. owstoni from L. c. canorum was always 1.0, regardless of the priors used in the analyses. In addition, the topology of the species tree inferred by BP&P was also consistent (94.2%-99.7%) with the one inferred by *BEAST. Therefore, our results suggested that L. c. canorum and L. c. owstoni could be regarded as two diverged evolutionary lineages.

17

3.3. Asymmetric post-divergence gene flow between L. c. canorum and L. c. owstoni The IMa analyses based on the 12 nuclear loci with and without CYTB revealed similar demographic patterns of the Hwamei complex (Table A.5), so we only focus on the results using nuclear DNA. The gene flow since the split of L. c. canorum and L. c. owstoni was asymmetrical, with the migration rate (per 1000 generations per gene) from L. c. canorum to L. c. owstoni (m2) around 0.0017 (95% CI = 0.0004-0.0076), whereas zero vice versa (m1= 0; 95% CI = 0.0001-0.0083; Fig. 4, Table A.5). The divergence time between L. c. canorum and L. c. owstoni was very recent (~ 210,000 years ago, 95% CI = 129,253-2,828,458), approximately the penultimate interglacial period and long before the last glacial maximum (LGM, ~26.5 Kya, Clark et al., 2009). The estimated effective population size (Ne) of L. c. canorum was 268,354 (95% CI=195,220-388,670), six times larger than that of L. c. owstoni (44,627, 95% CI = 24,575-79,228). The Ne of their most recent common ancestor was 98,495, with a larger uncertainty of estimation (95% CI=12,386-794,445, Fig. 4, Table A.5).

3.4. Different morphometrics between male L. c. canorum and L. c. owstoni Zheng et al. (1987) described L. c. owstoni as smaller than L. c. canorum, which is strongly supported by our data. All variables (except for BeL) of male Hwamei from Hainan were significantly shorter than those from the mainland (Table A.6). Results of discriminant analysis indicated that the five characteristics (BeL, WL, TL, TaL and aEL) can well distinguish specimens from Hainan and those from the mainland (Fig. A.2, the standardized coefficients of the linear combination of the five characteristics are listed in

18

Table A.7), despite the misclassification of several specimens [14.04% misclassification of specimens between Hainan and mainland, with most cases (five miss-assignments) occurring from Guangxi to Hainan, not shown].

4. Discussion Multiple pieces of genetic evidence show the divergence between L. c. owstoni and L. c. canorum. Our genetic data suggest that the split between the two Hwameis was fairly recent, approximately two glacial cycles ago. The speciation process between the two incipient Hwamei species did not fit a strictly allopatric speciation model: there was asymmetric post-divergence gene flow, probably through secondary contact during the last two glacial cycles when the sea level was low. Most of our genetic results (e.g., species tree, BP&P, multi-locus network, Structure) provide support to delimit L. c. owstoni from L. canorum. Moreover, morphometric analysis suggested a smaller body size and a shorter eyebrow length of male L. c. owstoni than those of L. canorum, despite the detection of historical gene flow since their split.

4.1. Species delimitation in Hwamei Evaluating the taxonomic status of species is critical to studies in macroecology (Meiri and Mace, 2007), phylogenetics (Dong et al., 2014), speciation (Li et al., 2010a), conservation (Isaac et al., 2004) and biogeography (Wiens, 2007). Even so, recognition of species is still challenging because of disagreement about the optimal criteria for defining them (e.g., de Queiroz, 1998; Harrison, 1998; Mayr, 1942). Previously, a species has been argued to equate to a metapopulation lineage that diverged from all other such

19

lineages (de Queiroz, 1998); the accumulation of concordant patterns among multiple independent characteristics (although not required, de Queiroz, 2007) has been thought to be robust evidence to distinguish such lineages, i.e. species delimitation (Avise and Ball, 1990; de Queiroz, 2007). However, recently diverged lineages pose the challenge of finding concordant evidence (Avise, 2000; de Queiroz, 2007). For example, a recent study of the light-vented/Taiwan bulbul (Pycnonotus sinensis/Pycnonotus taivanus) obtained conflicting results from various independent analyses based on different characteristics (Mckay et al., 2013). Both recent population isolation and continuous gene flow could lead to low levels of genetic divergence that would make species delimitation difficult (Petit and Excoffier, 2009). Based on a series of fifteen skins of L. c. owstoni from Hainan, Rothschild (1903) first described L. c. owstoni as a subspecies having a considerably paler plumage pattern than L. c. canorum from the adjacent continent. This could be partly supported by comparing images of Hwamei specimens taken from NZMC, IOZ, CAS (Fig. A.3). In the current analyses, the species tree of nuDNA that could control the effect of deep coalescence on gene tree topology supports the monophyly of L. c. ostwoni, although it is not consistent with the topology of the nuDNA gene tree based on the concatenated dataset. In addition, the multi-locus network, the Structure analysis and the coalescent species delimitation (*BEAST and BP&P) with nuclear loci all suggested the genetic distinctiveness of L.c. owstoni. Moreover, morphometric analyses (Table A.6) also indicated smaller body size and shorter eyebrow length of male L. c. owstoni than L. c. canorum from the adjacent mainland. According to the unified species concept (that species are segments of separately evolving metapopoulation lineages, even though the

20

operational criteria for delimiting lineages may differ) proposed by de Queiroz (2007), the taxonomic status of L. c. owstoni should be raised to full species rank, L. owstoni, even though its vocalization (see discussion below) may not be significantly different from that of L. canorum. However, the absence of any one or more properties (e.g., diagnosability, reproductive isolation) may not contradict a hypothesis of lineage separation (de Queiroz, 2007).

4.2. Incipient speciation with gene flow on the continental island According to the IMa estimates, L. canorum and L. owstoni began to diverge very recently, in the late Pleistocene (~0.21 Mya, million years ago). This is much more recent than a previous estimate using the CYTB gene alone (~0.6 Mya, Li et al., 2006). The potential overestimation of divergence dates by use of mitochondrial data alone has been documented in other studies (e.g., hu et al.,

1

Storchov et al., 2010). It seems that

such discrepancies might be partly caused by the fast coalescence rate of mitochondrial DNA, allowing the two monophyletic Hwamei groups to possess a most recent common mtDNA ancestor predating the date of their divergence. Moreover, since glaciations occurred every 100,000 years for the last one million years (reviewed in Williams et al., 1998), our genetic data suggest that L. owstoni split from its sister taxa only two glacial cycles ago. This indicates that Hainan, and maybe many other continental islands too, not only archives species originated millions of years ago, but also facilitates the genesis of new species. The IMa analyses revealed unambiguous signals of post-divergence gene flow between the mainland and island populations of Hwameis (Fig. 4). This suggests that the ancestor L. canorum was able to migrate southward and disperse into Hainan through 21

land bridges during glaciation, allowing secondary contact with L. owstoni and introgression. However, the low migration rate indicated that even though reproductive isolation was incomplete, the level of gene flow between two Hwameis was still restricted. Limited opportunity for secondary contact could have constrained the gene flow between two Hwamei lineages. However, Voris (2000) suggested that the north part of Hainan was connected with the nearby Asian continent more than half of the time for the last 250 thousand years, providing abundant opportunities for secondary contact. If true, alternative hypotheses such as directional selection (e.g., sexual selection, Gupta and Sundaran, 1994; Sawyer and Hartl, 1981) and other kinds of reproductive barriers (e.g., habitat barriers, Yu et al., 2000) might have been required to constrain gene flow between them. It has been suggested that male plumage and songs are major cues for species recognition and female choice in birds (reviewed by Price, 1998; Seddon et al., 2008). Although no comparison of L. canorum and L. owstoni song is available, individuals of L. owstoni always react to the playback of L. canorum when collecting samples. Such reaction suggests that L. owstoni can recognize vocal signals produced by L. canorum, so their vocalization might not be an effective cue for assortative mating between them. However, given the complexity of Hwamei’s songs (Tu and Severinghaus, 2004), it also possible that the preference of songs cannot be easily inferred through simple playback. Although tests of song reaction fall out of the scope of the current study, such analyses in the future might facilitate the understanding of their divergence process. The body size, eyebrow length and slight plumage color differences (Table A.6, Fig. A.3) between L. canorum and L. owstoni suggest that they might have contributed to sexual selection in

22

the Hwamei to some extent (although the effect of genetic drift on morphological difference cannot be ruled out). There are several potential causes for the asymmetric gene flow between the two Hwamei lineages, including 1) the larger Ne of L. canorum than L. owstoni, 2) specific reproductive mechanisms (such as asymmetric female preference towards L. canorum possibly because of its larger body size and darker plumage) that allowed asymmetric gene flow and 3) local selection (e.g., fit for warmer climate) for L. owstoni that might have inhibited gene flow from the island to the continent. Adaptation to different ecological or micro-ecological niches can create reproductive barriers between populations due to ecologically based divergent selection (Nosil, 2012). These could have influenced the level and direction of gene flow (e.g., Behm et al., 2010). The finding that L. owstoni may have speciated with gene flow has significant implications for the source of biodiversity of tropical islands (including the Sundaland, the world’s largest continental island system and one of the hotspots for biodiversity at risk, Myers et al., 2000). A recent phylogeographic study of the exceptionally rich mountain biota of Mt. Kinabalu in Sabah, East Malaysia, unambiguously demonstrates that a high proportion of endemic biota could have originated from regions outside Borneo (the eccentric endemics, Merckx et al., 2015). This is consistent with the expectation of an intermittent connection between Borneo, the adjacent islands and the continent in the Pleistocene ice ages (e.g., Voris, 2000). Speciation with gene flow for L. owstoni (this study) and the other lowland continental island endemics in Taiwan, another large continental island of East Asia (Dong et al., 2014; Hung et al., 2014; Li et al., 2010a), suggest that vicariance might not be essential for speciation of these eccentric

23

endemics; instead, non-allopatric speciation driven by various factors (e.g., selection, assortative mating, or even drift) might promote the biota genesis in these continental islands.

4.3. Incongruent topologies between the mitochondrial and nuclear gene trees The phylogenetic patterns of the mitochondrial and nuclear gene trees in the current study were incongruent, as was also found in previous studies (e.g., Hung et al., 2014; Mckay et al., 2013; Toews and Brelsford, 2012). In our study, the CYTB gene tree supported divergence between L. canorum and L. owstoni, while the concatenated nuclear gene tree provided an intermixed group of these two taxa. Such conflict is expected in the trees of recently diverged species (e.g., Harrington and Near, 2012; Hung et al., 2014; Mckay et al., 2013), and might be attributed to either deep coalescence (the coalescent time to the most recent common ancestor being longer than the history of a species; Maddison, 1997; Mckay et al., 2013) or introgressive hybridization (Harrington and Near, 2012). In theory, it would take in average four Ne generations for the nuclear DNA of a lineage to coalesce into a common ancestor back in time and reach monophyly (i.e., complete lineage sorting; Felsenstein, 2009; Kingman, 1982). However, the results of the IMa analysis indicated that the divergence time between L. owstoni and L. canorum was only around 210,000 years ago. Because the Ne of L. owstoni and L. canorum are estimated to be around 44 and 270 thousands individuals respectively, assuming a generation time of 2.5 years, the coalescent time of nuclear DNA would be more than 440 and 2,700 thousand years for the L. owstoni and L. canorum lineages. This is far

24

longer than the estimated divergence time between the two Hwameis. Such deep coalescent time could result in non-monophyly of the gene tree (e.g., Harrington and Near, 2012; Maddison, 1997) as we observed in the ML tree of nuclear genes. However, the Ne of mtDNA is only one quarter of that of nuclear DNA (Felsenstein, 2009). Therefore, the species’ mtDNA tree would be more likely to reach complete lineage sorting (e.g., Harrington and Near, 2012). The effect of recent divergence on the topology of the unresolved nuclear gene tree is supported by the results of the species tree (*BEAST and BP&P) analysis: when taking the effect of population demography (divergence time and Ne) into account, both L. owstoni and L. canorum are exclusive monophyletic clades, and are each other’s sister groups. In addition, introgression has been referred to as one of the potential weakness of relying largely on mitochondrial loci for species delimitation analyses (Bossu and Near, 2009, Harrington and Near, 2012); it may lead to a perceived genetic similarity even among phenotypically divergent species. With the observed monophyly in the mitochondrial gene tree of the three Hwamei taxa, no mitochondrial introgression was detected in our study (Fig. 2). The intermixed phylogenetic relationship between L. owstoni and L. canorum found by the analysis of concatenated nuclear genes cannot exclude the possibility of introgression in the nuclear DNA (e.g., possibly through sexual selection for male L. canorum from mainland). Even so, the coalescent analyses with multiple nuclear loci provide substantial support for recognition of three phylogenetic species in the Hwamei complex, outweighing the possibility of introgression through human interference.

25

4.4. Conservation of L. owstoni; contribution to the biodiversity of tropical and subtropical Asia It has been argued that species richness has been largely underestimated in tropical and subtropical Asia (Dong et al., 2014; Hung et al., 2014). Intensive human activities and animal trade in this region may also have influenced the process of species diversification [e.g., inter-specific gene flow can affect species integrity (Petit and Excoffier, 2009)], which is particularly severe for lineages that only diverged recently and have not arrived at complete reproductive isolation (e.g., Tu and Severinghaus, 2004). Hwamei have long suffered from trading. Large numbers of L. canorum, mainly from the Guangdong and Guangxi provinces of China, were caught and imported to Hainan for entertainment (personal communication with bird traders). Although no hybrids were found in the Structure analyses in the current study, given the high percentage (20.3%) of hybrids among L. taewanum (Li et al., 2010b), the recent divergence and cross vocal reaction of L. canorum and L. owstoni make it is reasonable to suspect an even higher potential for hybridization between them. Since the genetic diversity of L. owstoni is lower than that of the other two Hwamei groups, its wild population should be protected to avoid the loss of the unique local genetic material to hybridization (Moritz, 1994, reviewed by Rhymer and Simberloff, 1996). Legally regulated trading might be used to minimize the genetic admixture of L. owstoni with the exotic L. canorum.

Acknowledgements

26

We would like to thank Feng Dong and Chih-Ming Hung for advice on network and IMa analyses; Xiaolin Liao provided great help with the multivariate analyses. Special thanks are given to Fumin Lei, Qishan Wang, Ping Ding, Xiaojun Yang, Haiqing Tan and Zuohuai Wang who helped collect most of the Hwamei samples in the East Asian mainland and Hainan Island. Two anonymous reviewers, Scott Edwards and members of the labs of Rebecca T. Kimball and Edward L. Braun provided useful comments on revising the paper. Shuru Long and Danqing Zhang helped with some PCR experiments. We are grateful to the collaboration of the National Zoological Museum of China, Institute of Zoology, Chinese Academy of Sciences, the Museum of Biology, Sun yat-sen University, the Jianfengling National Nature Reserve and the Diaoluoshan National Nature Reserve. Chenxi Jia helped organize the measurement of specimens. We are also grateful to Alan Watson who had greatly improved the readability of this manuscript. This study was funded by the National Natural Science Foundation of China (grants 31360510 to N.W., 31301894 to B.L., 31260518 to J.W., and 31272328 and 31472013 to W.L.) and the Specimen Platform of China (2005DKA21403-JK).

Appendix A. Supplementary tables A.1-A.7 and figures A.1-A.3. Appendix B. Supplementary tree file.

References

Avise, J.C., 2000. Phylogeography: The History and Formation of Species. Harvard University Press, Cambridge, Massachusetts.

27

Avise, J.C., Ball, R., 1990. Principles of genealogical concordance in species concepts and biological taxonomy. Oxf. Surv. Evol. Biol. 7, 45–67. Barraclough, T.G., Vogler, A.P., 2000. Detecting the geographical pattern of speciation from species-level phylogenies. Am. Nat. 155, 419–434. Barrett, C.F., Freudenstein, J.V., 2011. An integrative approach to delimiting species in a rare but widespread mycoheterotrophic orchid. Mol. Ecol. 20, 2771-2786. Bank, C., Burger, R., Hermisson, J., 2012. The limits to parapatric speciation: dobzhanskymuller incompatibilities in a continent-island model. Genetics 191, 845-863. Behm, J.E., Ives, A.R., Boughman, J.W., 2010. Breakdown in postmating isolation and the collapse of a species pair through hybridization. Am. Nat. 175, 11–26. Bossu, C.M., Near, T.J., 2009. Gene trees reveal repeated instances of mitochondrial DNA introgression in orangethroat darters (Percidae: Etheostoma). Syst. Biol. 58, 114–129. Bruen, T.C., Philippe, H., Bryant, D., 2006. A simple and robust statistical test for detecting the presence of recombination. Genetics 172, 2665–2681. Bryant, D., Moulton, V., 2004. Neighbor-Net: an agglomerative method for the construction of phylogenetic networks. Mol. Biol. Evol. 21, 255–265. Campbell, N., Reece, J., 2002. Biology (6th ed.). Benjamin Cummings. Chu, J.H., Wegmann, D., Yeh, C.F., Lin, R.C., Yang, X.J., Lei, F.M., Yao, C.T., Zou, F.S., Li, S.H., 2013. Inferring the Geographic Mode of Speciation by Contrasting Autosomal and Sex-Linked Genetic Diversity. Mol. Biol. Evol. 30, 2519–2530. Clark, P.U., Dyke, A.S., Shakun, J.D., Carlson, A.E., Clark, J., Wohlfarth, B., Mitrovica, J.X., Hostetler, S.W., McCabe, A.M., 2009. The Last Glacial Maximum. Science 325, 710–714.

28

Clements, J.F., Schulenberg, T.S., Iliff, M.J., Roberson, D., Fredericks, T.A., Sullivan, B.L., Wood, C.L., 2015. The eBird/Clements checklist of birds of the world: v2015. Downloaded from http://www.birds.cornell.edu/clementschecklist/download/. Comes, H.P., Tribsch, A., Bittkau, C., 2008. Plant speciation in continental island floras as exemplified by Nigella in the Aegean Archipelago. Phil. Trans. R. Soc. B 363, 3083– 3096. Corstanje, R., Portier, K.M., Reddy, K.R., 2009. Discriminant analysis of biogeochemical indicators of nutrient enrichment in a Florida wetland. Eur. J. Soil Sci. 60, 974–981. Coyne, J.A., Orr, H.A., 2004. Speciation. Sinauer Associates, Sunderland, Massachusetts. Darwin, C., 1859. On the Origin of Species. John Murray, London. de Queiroz, K., 1998. The general lineage concept of species, species criteria, and the process of speciation: a conceptual unification and terminological recommendations, in: Howard, D.J., Berlocher, S.H. (Eds.), Endless Forms: Species and Speciation, Oxford University Press, New York, pp. 57–75. de Queiroz, K., 2007. Species concepts and species delimitation. Syst. Biol. 56, 879–886. Dong, F., Li, S.H., Zou, F.S., Lei, F.M., Liang, W., Yang, J.X., Yang, X.J., 2014. Molecular systematics and plumage coloration evolution of an enigmatic babbler (Pomatorhinus ruficollis) in East Asia. Mol. Phylogenet. Evol. 70, 76–83. Drummond, A.J., Rambaut, A., 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7, 214. Drummond, A.J., Ho, S.Y., Phillips, M.J., Rambaut, A., 2006. Relaxed phylogenetics and dating with confidence. PLoS Biol. 4, e88.

29

Earl, D.A., vonHoldt, B.M., 2012. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genet. Resour. 4, 359–361. Evanno, G., Regnaut, S., Goudet, J., 2005. Detecting the number of clusters of individuals using the software structure: a simulation study. Mol. Ecol. 14, 2611–2620. Excoffier, L., Smouse, P.E., Quattro, J.M., 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131, 479–491. Excoffier, L., Laval, G., Schneider, S., 2005. ARLEQUIN (version 3.01): an integrated software package for population genetic data analysis. Evol. Bioinform. 1, 47–50. Falush, D., Stephens, M., Pritchard, J.K., 2003. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164, 1567–1587. Felsenstein, J., 2009. Theoretical Evolutionary Genetics, pp 364-368. Seatle, Washington. Garant, D., Kruuk, L.E.B., Wilkin, T.A., McCleery, R.H., Sheldon, B.C., 2005. Evolution driven by differential dispersal within a wild bird population. Nature 433, 60–65. Gaston, K.J., 2003. The Structure and Dynamics of Geographic Ranges. Oxford University Press, Oxford. Gavrilets, S., Vose, A., 2005. Dynamic patterns of adaptive radiation. Proc. Natl Acad. Sci. USA 102, 18040–18045. Gupta, J.P., Sundaran, A.K., 1994. Some evidence of incipient speciation in Drosophila kikkawai. Genome 37, 1041–1044.

30

Hackett, S.J., Kimball, R.T., Reddy, S., Bowie, R.C.K., Braun, E.L., Braun, M.J., Chojnowski, J.L., Cox, W.A., Han, K.-L., Harshman, J., Huddleston, C.J., Marks, B.D., Miglia, K.J., Moore, W.S., Sheldon, F.H., Steadman, D.W., Witt, C.C., Yuri, T., 2008. A phylogenomic study of birds reveals their evolutionary history. Science 320, 1763–1768. Harrigan, R.J., Mazza, M.E., Sorrenson, M.D., 2008. Computation vs. cloning: evaluation of two methods for haplotype determination. Mol. Ecol. Resour. 8, 1239–1248. Harrison, R.G., 1998. Linking evolutionary pattern and process, in: Howard, D.J., Berlocher, S.H. (Eds.), Endless Forms: Species and Speciation, Oxford University Press, New York, pp.19–31. Harrington, R.C., Near, T.J., 2012. Phylogenetic and coalescent strategies of species delimitation in snubnose darters (Percidae: Etheostoma). Syst. Biol. 61, 63-79. Heled, J., Drummond, A.J., 2010 Bayesian inference of species trees from multilocus data. Mol. Biol. Evol. 27, 570–580. Hey, J., Nielsen, R., 2007. Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proc. Natl Acad. Sci. USA 104, 2785–2790. Hillier, L.W., Miller, W., Birney, E., Warren, W., et al., 2004. Sequence and comparative analysis of the chicken genome provide unique perspectives on vertebrate evolution. Nature 432, 695–716. Hoskin, C.J., Higgie, M., McDonald, K.R., Moritz, C., 2005. Reinforcement drives rapid allopatric speciation. Nature 437, 1353–1356. Howard, R., Moore, A., 1994. A Complete Checklist of the Birds of the World, 2nd edn. Academic Press, San Diego, CA.

31

Hudson, R.R., Kreitman, M., Aguade, M., 1987. A test of neutral molecular evolution based on nucleotide data. Genetics 116, 153–159. Hudson, R.R., Slatkin, M., Maddison, W.P., 1992. Estimation of levels of gene flow from DNA-sequence data. Genetics 132, 583–589. Hung, C.M., Hung, H.Y., Yeh, C.F., Fu, Y.Q., Chen, D., Lei, F., Yao, C.T., Yao, C.J., Yang, X.J., Lai, Y.T., Li, S.H., 2014. Species delimitation in the Chinese bamboo partridge Bambusicola thoracica (Phasianidae; Aves). Zool. Scr. 43, 562–575. Huson, D.H., Bryant, D., 2006. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23, 254–267. Isaac, N.J.B., Mallet, J., Mace, G.M., 2004, Taxonomic inflation: its influence on macroecology and conservation. Trends Ecol. Evol. 19, 464–469. Joly, S., Bruneau, A., 2006. Incorporating allelic variation for reconstructing the evolutionary history of organism from multiple genes: an example from Rosa in North America. Syst. Biol. 55, 623–636. Katoh, K., Asimenos, G., Toh, H., 2009. Multiple alignment of DNA sequences with MAFFT. Methods Mol. Biol. 537, 39–64. Kingman, J.F.C., 1982. The coalescent. Stochastic Process. Appl. 13, 235–248. Lambeck, K., Chappell, J., 2001. Sea level change through the last glacial cycle. Science 292, 679–686. Lande, R., 1980. Genetic Variation and Phenotypic Evolution During Allopatric Speciation. Am. Nat. 116, 463–479. Leaché, A.D., Fujita, M.K., 2010. Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus). P. Roy. Soc. Lond. B Biol. Sci. 277, 3071–3077.

32

Leaché, A.D., Harris, R.B., Rannala, B., Yang, Z., 2014. The influence of gene flow on species tree estimation: a simulation study. Syst. Biol. 63, 17–30. Li, S.H., Li, J.W., Han, L.X., Yao, C.T., Shi, H., Lei, F.M., Yen, C., 2006. Species delimitation in the Hwamei Garrulax canorus. Ibis 148, 698–706. Li, J.W., Yeung, C.K.L., Tsai, P.W., Lin, R.C., Yeh, C.F., Yao, C.T., Han, L., Hung, L.M., Ding, P., Wang, Q., Li, S.H., 2010a. Rejecting strictly allopatric speciation on a continental island: prolonged postdivergence gene flow between Taiwan (Leucodioptron taewanus, Passeriformes Timaliidae) and Chinese (L. canorum canorum) hwameis. Mol. Ecol. 19, 494–507. Li, S.H., Yeung, C.K.L., Han, L., Le, M.H., Wang, C.X., Ding, P., Yao, C.T., 2010b. Genetic introgression between an introduced babbler, the Chinese hwamei Leucodioptron c. canorum, and the endemic Taiwan hwamei L. taewanus: a multiple marker systems analysis. J. Avian Biol. 41, 64–73. Liao, X., Inglett, P.W., Inglett, K.S., 2014. Vegetation and microbial indicators of nutrient status: Testing their consistency and sufficiency in restored calcareous wetlands. Ecol. Indic. 46, 358–366. Librado, P., Rozas, J., 2009. DnaSP v5: a software for compre- hensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. Lischer, H.E.L., Excoffier, L., 2012. PGDSpider: An automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics 28, 298-299. MacKinnon, J., Phillipps, K., 2000. A Field Guide to the Birds of China. Oxford University Press. Maddison, W.P., 1997 Gene trees in species trees. Syst. Biol. 46, 523–536.

33

Mayr, E., 1942. Systematics and the Origin of Species from the Viewpoint of a Zoologist. Columbia University Press, New York. Merckx, V.S., Hendriks, K.P., Beentjes, K.K., Mennes, C.B., Becking, L.E., Peijnenburg, K.T., Afendy, A., Arumugam, N., de Boer, H., Biun, A., Buang, M.M., 2015. Evolution of endemism on a young tropical mountain. Nature 524, 347–350. McKay, B.D., Mays Jr, H.L., Wu, Y., Li, H., Yao, C-T., Nishium, I., Zou, F., 2013. An empirical comparison of character-based and coalescent-based approaches to species delimitation in a young avian complex. Mol. Ecol. 22, 4943–4957. Meiri, S., Mace, G.M., 2007. New taxonomy and the origin of species. PLoS Biol. 5, e194. Moritz, ., 1994. Defining ‘evolutionary significant units’ for conservation. Trends. Ecol. Evol. 9, 373–375. Myers, N., Mittermeier, R.A., Mittermeier, C.G., da Fonseca, G.A.B., Kent, J., 2000. Biodiversity hotspots for conservation priorities. Nature 403, 853–858. Nei, M., 1987. Molecular Evolutionary Genetics. New York, U.S.: Columbia Univ. Press. Niemiller, M., Fitzpatrick, B., Miller, B., 2008. Recent divergence with gene flow in Tennessee cave salamanders (plethodontidae: gyrinophilus) inferred from gene geneaologies. Mol. Ecol. 17, 2258–2275. Nosil, P., 2012. Ecological speciation. OUP Oxford. Petit, R.J., Excoffier, L., 2009. Gene flow and species delimitation. Trends Ecol. Evol. 24, 386–393. Posada, D., Crandall, K.A., 1998. MODELTEST: testing the model of DNA substitution. Bioinformatics. 14, 817–818.

34

Price, T., 1998. Sexual selection and natural selection in bird speciation. Philos. Trans. R. Soc. Lond. B Biol. Sci. 353, 251–260. Pritchard, J.K., Stephens, M., Donnelly, P., 2000. Inference of population structure using multilocus genotype data. Genetics 155, 945–959. Rambaut, A., Drummond, A.J., 2007. Tracer v1.4. Available at: http://tree.bio.ed.ac.uk/software/tracer/. Rhymer, J.M., Simberloff, D., 1996. Extinction by hybridization and introgression. Annu. Rev. Ecol. Syst. 27, 83–109. Ronquist, F., Huelsenbeck, J.P., 2003. MrBayes 3: Bayesian phy- logenetic inference under mixed models. Bioinformatics 19, 1572–1574. Rosenberg, N.A., 2004. DISTRUCT: a program for the graphical display of population structure. Mol. Ecol. Notes 4, 137–138. Rothschild, H.W., 1903. [Mr. Rothschild also made remarks on a large collection of birds received from the Island of Hainan.] Bull. Brit. Orn. Club 14, 6–9. Sawyer, S., Hartl, D., 1981. On the evolution of behavioral reproductive isolation: The Wallace effect. Theor. Popul. Biol. 19, 261–273. Seddon, N., Merrill, R.M., Tobias, J.A., 2008. Sexually selected traits predict patterns of species richness in a diverse clade of suboscine birds. Am. Nat. 171, 620–631. Shaner, P.J.L., Tsao, T.H., Lin, R.C., Liang, W., Yeh, C.F., Yang, X.J., Lei, F.M., Zhou, F., Yang, C.C., Hung, L.M., Hsu, Y.C., Li, S.H., 2015. Climate niche differentiation between two passerines despite ongoing gene flow. J. Anim. Ecol. 84, 829–839. Shi, M., Chen, C., Xu, Q., Lin, H., Liu, G., Wang, H., Wang, F., Yan, J., 2002. The role of Qiongzhou Strait in the seasonal variation of the South China Sea circulation. J. Phys. Oceanogr. 32, 103-121. 35

Slatkin, M., 1987. Gene flow and the geographic structure of natural populations. Science 236, 787–792. Sorenson, M.D., Fleischer, R.C., 1996. Multiple independent transpositions of mitochondrial DNA control region sequences to the nucleus. Proc. Natl Acad. Sci. USA 93, 15239– 15243. Stamatakis, A., 2006. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics 22, 2688– 2690. Stephens, M., Smith, N.J., Donnelly, P., 2001. A new statistical method for haplotype reconstruction from population data. Am. J Hum. Genet. 68, 978–989. Storchová, R., Reif, J., Nachman, M.W., 2010. Femaleheterogametyand speciation: reduced introgression of the Z chromosome between two species of nightingales. Evolution 64, 456–471. Swinhoe, R., 1859. Notes on some new species of birds found on the Island of Formosa. J. N. China Br. Roy. Ass. Soc. 1, 228. Swofford, D.L., 2002. PAUP*, Phylogenetic analysis using parsimony (* and other methods). Version 4. Sinauer Associates, Sunderland, Massachusetts. Tajima, F., 1989. Statistical method for testing the neutral mutation hypothesis by 913. DNA polymorphism. Genetics 123, 585–595. Toews, D.P.L., Brelsford, A., 2012. The biogeography of mitochondrial and nuclear discordance in animals. Mol. Ecol. 21, 3907–3930. Tu, H.W., Severinghaus, L.L., 2004. Geographic variation of the highly complex Hwamei (Garrulax canorus) songs. Zool. Stud. 43, 629–640.

36

Voris, H.K., 2000. Maps of Pleistocene sea levels in Southeast Asia: shorelines, river systems and time durations. J Biogeogr. 27, 1153–1167. Warren, W.C., Clayton, D.F., Ellegren, H., Arnold, A.P., et al., 2010. The genome of a songbird. Nature 464, 757–762. Weir, J.T., Schluter, D., 2008. Calibrating the avian molecular clock. Mol. Ecol. 17, 2321– 2328. Weir, B.S., Cockerham, C.C., 1984. Estimating F-statistics for the analysis of population structure. Evolution 38, 1358–1370. Whittaker, . ., ern ndez-Palacios, J.M., 2007. Island biogeography: ecology, evolution, and conservation, 2nd edn. Oxford, UK: Oxford University Press. Wiens, J.J., 2007. Species delimitation: new approaches for discovering diversity. Syst. Biol. 56, 875–878. Williams, M.A.J., Dunkerley, D.L., de Deckker, P., Kershaw, A.P., Chappel, J., 1998. Quaternary Environments. Arnold, London. Won, Y.J., Hey, J., 2005. Divergence population genetics of chimpanzees. Mol. Biol. Evol. 22, 297–307. Xing, F.W., Wu, T.L., Li, Z.X., Ye, H.G., Chen, B.H., 1995. Endemic plants of Hainan island. J Trop. Subtrop. Bot. 3, 1–12. Yan, J., 2006. Paleontology and ecologic environmental evolution of the Quaternary in Hainan Island. J Palaeogeogr. 8, 103–115. Yang, Z., Rannala, B., 2010. Bayesian species delimitation using multilocus sequence data. Proc. Natl. Acad. Sci. USA 107, 9264–9269. Yang, Z., Rannala, B., 2014. Unguided species delimitation using DNA sequence data from

37

multiple loci. Mol. Biol. Evol. 31, 3125–3135. Yang, Z., 2015. The BPP program for species tree estimation and species delimitation. Curr. Zool. 61, 854–865. Yeung, C.K.L., Tsai, P.W., Chesser, R.T., Lin, R.C., Yao, C.T., Tian, X.H., Li, S.H., 2011. Testing founder effect speciation: divergence population genetics of the spoonbills Platalea regia and Pl. minor (Threskiornithidae, Aves). Mol. Biol. Evol. 28, 473–482. Yu, G., Chen, X., Ni, J., Cheddadi, R., et al. 2000. Palaeovegetation of China: a pollen databased synthesis for the mid-Holocene and last glacial maximum. J. Biogeogr. 27, 635– 664. Zheng, G.M., 2011. A checklist on the classifircation and distribution of the birds of China (second edition). Science Press, Beijing. Zheng, Z., Long, Z., Zheng, B., 1987. Aves, Vol. 11, Fauna Sinica, Science Press, Beijing.

38

Table 1 Genetic characteristics of the mtDNA and 12 nuclear genes in Hwamei.

Locus/gene CYTB

ABCC5

ADIPOR1

ALDOB

DPYSL3

F1N9H4

HMGB2

HPGDS

LMNA

Populations Hainan Mainland Taiwan Hainan Mainland Taiwan Hainan Mainland Taiwan Hainan Mainland Taiwan Hainan Mainland Taiwan Hainan Mainland Taiwan Hainan Mainland Taiwan Hainan Mainland Taiwan Hainan

N1 14 29 12 12 67 15 13 32 26 13 79 16 13 73 25 13 74 20 13 77 18 12 69 10 13

Inheritance2 Mitochondrial

Length (bp) 1143

Chromosome 9

551

Chromosome 26

Chromosome Z

360

574

Chromosome 13

696

Chromosome 4

386

Chromosome 4

Chromosome 4

Chromosome 2

505

772

252

Alleles 6 23 7 6 17 4 2 5 1 4 1 1 1 14 3 2 18 1 3 9 2 5 23 1 4

S3 9 32 10 4 9 3 1 3 0 3 0 0 0 10 2 1 17 0 2 7 2 9 31 0 2

π4 × E-3 1.85 3.27 2.85 2.91 3.26 1.26 0.22 1.58 0.00 0.76 0.00 0.00 0.00 1.52 0.12 0.38 2.48 0.00 0.69 1.75 0.22 5.28 7.39 0.00 4.11

θw5 × E-3 2.48 7.13 2.90 1.94 2.99 1.37 0.75 1.82 NA 1.37 NA NA NA 2.59 0.64 0.68 7.91 NA 1.04 2.47 0.96 3.13 7.56 NA 2.08

Dxy6 0.0149

Phi7 Model8 0.5788 TVM+I+G

0.0319 0.0042

0.1207 K80+I

0.0031 0.0012

0.8633 TrN

0.0069 0.0004

1

TrNef

0.0000 0.0009

0.8386 HKY+I

0.0033 0.0015

0.2713 K81uf+I

0.0058 0.0018

0.2225 TVM

0.0049 0.0096

0.2205 HKY+I

0.0074 0.0048

0.2767 HKY+I

39

Mainland 76 Taiwan 7 MOCS1 Hainan 11 Chromosome 3 590 Mainland 73 Taiwan 12 SDC2 Hainan 12 Chromosome 2 683 Mainland 79 Taiwan 6 TLE4 Hainan 13 Chromosome Z 691 Mainland 79 Taiwan 10 VLDLR9 Hainan 13 Chromosome Z 415 Mainland 78 Taiwan 15 The complete gene name can be found in Appendix A, Table A.2.

6 2 2 10 1 3 10 2 1 4 1 5 18 1

3 1 1 9 0 2 9 1 0 3 0 3 12 0

3.14 2.14 0.29 2.83 0.00 0.63 0.37 0.44 0.00 0.05 0.00 3.48 4.65 0.00

2.13 1.25 0.46 2.74 NA 0.78 2.34 0.48 NA 0.77 NA 1.89 5.14 NA

0.0028 0.0034

0.0704 F81

0.0060 0.0006

0.7454 HKY+I

0.0004 0.0000 0.0000 0.0045

1

HKY

0.4218 HKY

0.0079

1

N: the number of individuals sequenced; 2 The chromosome number refers to the chicken genome.

3

S: the number of segregating sites; 4π: nucleotide diversity; 5Watterson’s θ (θw): the genetic diversity; 6Dxy: genetic distance between

L. c. canorum and L. c. owstoni and between L. c. canorum and L. taewanum; 7

Phi: recombination test indicated no evidence of recombination within locus; 8Model: the optimal substitution model chosen by AIC

criterion.

40

Table 2 Estimation of the substitution rate at each nuclear locus relative to that of CYTB based on the average genetic distance between each Hwamei species/subspecies and the outgroup. Net genetic distance

1

Locus

owstonioutgroup

canorumoutgroup

taewanumoutgroup

Average net distance

CYTB ABCC5 ADIPOR1 ALDOB DPYSL3 F1N9H4 HMGB2 HPGDS LMNA MOCS1 SDC2 TLE4 VLDLR9

0.087 0.014 0.027 NA 0.022 0.021 0.017 NA 0.023 NA 0.004 0.009 0.037

0.084 0.015 0.027 NA 0.021 0.021 0.016 NA 0.023 NA 0.004 0.009 0.037

0.091 0.017 0.027 NA 0.021 0.027 0.018 NA 0.022 NA 0.004 0.009 0.030

0.087 0.015 0.027 NA 0.021 0.023 0.017 NA 0.023 NA 0.004 0.009 0.035

Relative to net distance at CYTB 1 0.172 0.310 NA 0.241 0.264 0.195 NA 0.264 NA 0.046 0.103 0.402

s/s/y1 ×10-8

s/l/y2 ×10-8

1.035 0.178 0.321 NA 0.250 0.274 0.202 NA 0.274 NA 0.048 0.107 0.416

1183.0 98.3 115.6 NA 173.9 105.6 102.1 NA 69.0 NA 32.5 74.0 172.8

lower ranges

higher ranges

s/s/y ×10-8

s/l/y ×10-8

s/s/y ×10-8

s/l/y ×10-8

0.600 0.103 0.186 NA 0.145 0.159 0.117 NA 0.159 NA 0.028 0.062 0.241

685.8 57.0 67.0 NA 100.8 61.2 59.2 NA 40.0 NA 18.9 42.9 100.2

1.505 0.259 0.467 NA 0.363 0.398 0.294 NA 0.398 NA 0.069 0.156 0.605

1720.2 143.0 168.1 NA 252.8 153.6 148.5 NA 100.3 NA 47.3 107.6 251.3

s/s/y indicates substitution/site/year; 2s/l/y indicates substitution/locus/year.

41

Table 3 FST comparison between Hwamei groups using either the CYTB gene or the concatenated nuclear loci. Population pairwise FST CYTB Hainan Mainland Taiwan Outgroup Nuclear Hainan Mainland Taiwan Outgroup

Hainan

Mainland Taiwan

0.81258 0.92682 0.97387

0.90128 0.95834 0.96621

0.86261 1 1

0.83543 0.98007

1

42

Figure Legends:

Figure 1. Sampling areas for different Hwamei populations. The approximate range of the Chinese Hwamei (dashed line) was drawn according to Birdlife International (http://www.birdlife.org/datazone/speciesfactsheet.php?id=32442). Abbreviations of sampled populations are as follows: SX, Shannxi Province; AH, Anhui Province; ZJ, Zhejiang Province; GZ, Guizhou Province; HN, Hainan Province; TW, Taiwan; VN, Viet Nam.

Figure 2. Partitioned maximum likelihood analyses of Hwamei based on A: the concatenated nuclear + CYTB gene and B: the concatenated nuclear genes. Bootstrap values of some focal clades are indicated on the nodes. Individuals colored in blue are from L. c. canorum.

Figure 3. A: Maximum Likelihood analysis of Hwamei in RAxML based on the single CYTB gene. Bootstrap supports (above branch at node) and posterior probability (below branch at node) are provided for the divergence of major lineages. B: Multi-locus genetic network based on 12 nuclear loci. C: Species tree of the Hwamei complex inferred from *BEAST analyses. Posterior probabilities of clade support are labeled at each node. D: Bayesian Structure clustering results based on nuclear sequences of Hwamei, indicating three genetic clusters.

Figure 4. IMa estimates of demographic parameters for the divergence history between L.

43

canorum and L. owstoni based on 12 nuclear loci. A: Effective population sizes (Ne); B: the migration rate (average number of migrations per 1000 generations per gene copy) from L. canorum to L. owstoni (m2) and vice versa (m1); and C: the divergence time in million years (T) between the two taxa.

44

Figure 1

SX AH ZJ GZ TW VN HN

Figure 2

76 100 32

C823 C581 C820 C408

T1141 T3443 T2782 T3120 T2677 T800 T3118 T2689 T1097 T3486 T774 T772 T3441 T3113 T773 T280 T3117 T2688 T2681 T2684 T2692 T3116 T3437 T3485 T277 T541 T1099 T3500 T3111 T770

C4364 C590 C342

100 L.taewanum

18

0.005 A: Partitioned ML (12 Nuclear loci+CYTB)

C821

C4363 C346 C365 C814 C815 C367 C587 C578 C411 C400 C575 C819 C502 C352 C4362 C387 C497 C787 C500 C498 C493 C330 C362 C577 C584 C822 C4086 C586 C392 C491 C495 C499 C591 C816 C358 C4360 C403 C490 C796 C379 C355 C576 C4466 C580 C395 C4148 C501 C350 C783 C504 C4361 C782 C496 C588 C492 C503 C592 C397 C4144 C393 C585 C39 C4359 C582 C370 C579 C817 H14 H6 L.c.owstoni H11 C820 C581 C408 C818 C823 T770 T2782 T2688 T1141 T2692 T3120 T3486 T277 T1097 T3117 T3116 T2681 T2689 T800 T772 T774 T3441 T280 T3113 T773 T3443 T2677 T3118 T2684 T3111 T3500 T3437 T3485 T1099 T541 C399 C489 H3 H2 H4 H10 H7 L.c.owstoni H12 H5 H8 H9 H1

71

25 0.004

B: Partitioned ML (12 Nuclear loci)

L.taewanum

C352 C393 C585 C39 C582 C370 C4359 C330 C387 C497 C502 C493 C4362 C403 C4360 C796 C362 C4148 C358 C395 C576 C577 C818 C495 C579 C4466 C816 C350 C501 C503 C592 C4361 C580 C787 C782 C490 C379 C492 C496 C588 C392 C822 C584 C586 C4144 C4086 C491 C499 C591 C500 C498 C504 C783 C575 C819 C815 C817 C4364 C4363 C821 C411 C578 C587 C400 C355 C342 C583 C590 C399 C365 C397 C814 C346C367 C489 H14 H6 H11 H2 H3 H12 L.c.owstoni H9 H8 H4 H10 H5 H7 H1

C583

Figure 3

96 1

A

100

C230 C248 C243 C374 C244 C39 C578 C414 C330 40 C353 0.58 C574 C416 C593 C569 C415 C417 C342 C341 C380 C242 C232 41 C231 C352 0.49 C373 C382 C381 C228 C371 C372 H951 H777 H12 96 H13 0.99 H779 H8 H14 H5 H949 H6 H2 H11 H9 H10 H950 H3 H4 H7 H1

T801 T773 T802 T280 T340 T279 T768 T277 T276 T799 T769 T604

L.taewanum

B

L.c.canorum

0.02 L.c.owstoni L.c.canorum L.taewanum Ga

r ru

L.c.owstoni

G.monileger

lax

mo ni leg er

G.monileger

0.02 1

L.taewanum

1

L.c.owstoni L.c.owstoni

C

D

1 0.004

L.c.canorum

L.c.canorum

L.taewanum

Figure 4

3.5

5

0.8

owstoni

3.0

0.7

2.5

0.6 0.5

2.0

m1: owstoni

canorum

m2: canorum

owstoni

4 3

0.4

1.5

ancestor canorum

1.0

2

0.3 0.2

0.5

1

0.1 0

0 0

0.3 Ne

0.6 Millions

0 0

0.005

0.01 m

0.015

0.02

0

1

2 T

3

4 Millions

Graphical Abstract

96 1

A

100

C230 C248 C243 C374 C244 C39 C578 C414 C330 40 C353 0.58 C574 C416 C593 C569 C415 C417 C342 C341 C380 C242 C232 41 C231 C352 0.49 C373 C382 C381 C228 C371 C372 H951 H777 H12 96 H13 0.99 H779 H8 H14 H5 H949 H6 H2 H11 H9 H10 H950 H3 H4 H7 H1

T801 T773 T802 T280 T340 T279 T768 T277 T276 T799 T769 T604

L.taewanum

B

L.c.canorum

0.02 L.c.owstoni L.c.canorum L.taewanum Ga

r ru

L.c.owstoni

G.monileger

lax

mo ni leg er

G.monileger

0.02 1

L.taewanum

1

L.c.owstoni L.c.owstoni

C

D

1 0.004

L.c.canorum

L.c.canorum

L.taewanum

Highlights  Genetic and morphological evidences suggest divergence between L. c. owstoni and L. c. canorum.  Hainan Hwamei could be regarded as a full species (L. owstoni).  IMa analysis suggests that the two Hwameis split ~0.2 million years ago.  Limited asymmetrical gene flow probably existed during the divergence of the two Hwameis.  A strictly allopatric process can’t explain speciation of a continental island endemic species.

45