A multilocus molecular phylogeny for the avian genus Liocichla (Passeriformes: Leiothrichidae: Liocichla)

Background

The babblers are a morphologically and ecologically diverse group of some 300 species and over 800 named taxa of predominantly tropical Old World passerines, with the exception of the New World Wrentit (Chamaea fasciata) [1]. As a prominent component of the forest biotas in tropical and subtropical Asia, babblers are a potential model system for investigating ecological, evolutionary and biogeographic hypotheses [2]. However, developing a robust phylogeny for this diverse group is a necessary prerequisite for integrating the group into broader evolutionary and biogeographic studies. Given the number of species as well as their morphological and ecological diversity and high degree of character convergence this goal is not without its challenges but recent work has made considerable progress in imposing some taxonomic rigor on this diverse group [3-14].

Babblers have traditionally been assigned a family status (Timaliidae) but that taxonomic status has never rested on a solid foundation [2, 15-17]. The Timaliidae have been called a taxonomic "dustbin" [1] or a "scrap basket" [18] and recent phylogenetic analyses highlight the taxonomic problems associated with this group (reviewed by Alström et al. [4]). Stachyris and Alcippe are polyphyletic [7, 11, 19], Pomatorhinus and Actinodura are paraphyletic [8, 9, 11], Pteruthius and Erpornis are actually vireos, Vireonidae [6, 13, 20], Pnoepyga has been suggested to be placed in a separate family, Pnoepygidae [4, 5, 10, 21] and the traditional babbler genus Yuhina is allied with the Zosteropidae [7, 10, 11]. Alström et al. [5] supports the removal of Elachura formosa (previously Spelaeornis formosus) not just from an affiliation with the babblers but from the Sylvioidea superfamily entirely and into its own monotypic family, Elachuridae. Gelang et al. [10] has recognized five clades encompassing the babblers with four assigned the rank of subfamily (Leiothrichinae, Pellorneinae, Timaliinae and Zosteropinae) within a monophyletic Timaliidae that is sister to Sylviidae. Fregin et al. [21], Alström et al. [4, 5] and Gill and Donsker [22] elevate each of the subfamilies within the Timaliidae sensu Gelang et al. [10] to the rank of family whereas Moyle et al. [11] treat Leiothrichinae, Pellorneinae, Timaliinae as subfamilies within Timaliidae and Zosteropidae and Sylviidae at the family level. While these studies have made great strides in elucidating the family-level relationships, the scope of the family has created a logistical hurdle to creating a comprehensive molecular phylogeny.

One babbler genus that has been incompletely sampled in prior studies is Liocichla. Five species are currently recognized in the genus: Liocichla phoenicea (Red or Crimson-faced Liocichla), L. ripponi (Scarlet-faced Liocichla), L. bugunorum (Bugun Liocichla), L. omeiensis (Grey-cheeked Liocichla) and L. steerii (Steere's Liocichla) [23, 24]. These species fall within the Leiothrichidae sensu (Gill and Donsker [22]) and are closely allied with the genera Actinodura, Heterophasia, Leiothrix, and Minla [6, 10, 11].

Three of the five species within the Liocichla are restricted range species (Fig. 1). L. bugunorum is a recently described taxon with a very limited range in the Northeast Indian state of Arunachal Pradesh [25]. L. omeiensis is a Chinese endemic also with a very limited range primarily confined to Sichuan Province [1, 24]. L. omeiensis is considered vulnerable with a decreasing population trend [26] and while data on L. bugunorum are currently insufficient for an accurate population assessment, its restricted range and limited number of sightings likely indicate that it too is vulnerable [1, 24, 27]. L. steerii is an island endemic confined to the island of Taiwan with plumage reminiscent of L. bugunorum and L. omeiensis [1, 24, 28].

Figure  1.  Map of ingroup sampling locations and distributions of taxa. Distributions based on Collar and Robson [1].

DownLoad: Full-Size Img PowerPoint

The more geographically widespread L. phoenicea has historically been split into four subspecies with L. p. phoenicea in the Eastern Himalayas and Northeast Assam, L. p. bakeri in Southern Assam, Western and Northern Myanmar and Northwestern Yunnan, L. p. ripponi in Eastern Myanmar, Northwestern Thailand and Southwestern Yunnan and L. p. wellsi in Southeastern Yunnan and Northern Indochina [29]. Collar's [23] examination of plumage characters and Dickinson and Christidis [24] elevates L. p. ripponi and L. p. wellsi to recognition as a single species (L. ripponi) and L. p. phoenicea and L. p. bakeri are likewise united into a single species (L. phoenicea). In addition to being distinguished by plumage L. phoenicea and L. ripponi sensu Collar [23] are also parapatric with L. phoenicea occupying the more Northwestern portion of the distribution of L. phoenicea sensu lato in Northern and Western Myanmar, Eastern India, Nepal and Western Yunnan and L. ripponi occupying South and Southeastern Myanmar, Southwestern Yunnan, Northwestern Thailand, Northeastern Guangxi and extreme North Vietnam and Laos with the two species' distributions approximately divided in Myanmar by the Myitnge River ([1], Fig. 1) albeit with potential exchange occurring across a narrow contact zone in this region of Eastern Myanmar [23].

No babbler molecular phylogeny to date has included all Liocichla species. Luo et al. [30] included L. steerii and L. omeiensis, Gelang et al. [10] and Dong et al. [9] each included only L. steerii, and Cibois [6] and Moyle et al. [11] included L. steerii and L. phoenicea (albeit in the case of L. phoenicea with samples from different populations that Collar [23] would delimit as belonging to different species). Generating a robust molecular phylogeny for Liocichla that includes all currently recognized taxa would be a valuable piece of the babbler puzzle and may provide a point of comparison for broader studies of speciation and biogeography in Asia. A reliable and comprehensive phylogeny for the babblers is particularly relevant for elucidating the history of continental island biotas in Asia such as those in Taiwan, Hainan, Indochina and the Greater Sunda islands and for the biotas associated with the uplift of the Himalayas [31-33].

Methods

Sampling and laboratory methods

We obtained DNA sequence data from all five species of Liocichla currently recognized by Collar and Robson [1] and Gill and Donsker [22]. Sequence data for cytb (KJ456321.1) and ND2 (KJ455478.1) derived from L. bugunorum were generated in a previous study [33] while the remainder of our samples originated from specimen vouchered frozen tissue samples. Samples from the genera Minla and Leiothrix, the putative closest extant relatives of Liocichla [11], were obtained and used as outgroups. Two samples of L. ripponi were derived from specimens of captive origin in the Louisiana State University Museum of Natural History (LSUMZ) collection originally determined as L. phoenicea (LSUB20702 and LSU37270). Only one of these samples was associated with a skin while the other was associated with skeletal material. The single skin available for the LSUMZ captive specimen exhibited diagnostic plumage consistent with L. ripponi (grey crown, scarlet face extending above the eye, grey underparts, terminal bar on the rectrices; [1, 23]. These two specimens from LSUMZ were nearly identical genetically to each other and nearly genetically identical to a field collected L. ripponi collected from China (KU10078). For two ingroup samples (L. phoenicea KU15195 and L. ripponi KU10078) and three outgroup samples (Leiothrix argentauris KU15203, Leiothrix lutea KU6724 and Minla ignotincta KU11346) sequence data were available from prior studies. We re-sequenced these same samples for this study and use these data in our analysis but cite the Genbank accession numbers from the prior work as our sequences were shorter and did not differ significantly from the original Genbank submissions (Additional file 1: Table S1). Voucher information for all samples is available in Table 1.

Table  1.  Associated data for tissue/blood samples for ingroup (Liocichla sp.) and outgroup (Leiothrix argentauris, Leiothrix lutea and Minla ignotincta) genetic samples

Catalog number	Loaning institution	Taxon	Country	Locality	Latitude	Longitude
DT30 (feather sample)	Wildlife Institute of India, Dehradun	Liocichla bugunorum	India	Eaglenest Wildlife Sanctuary, Arunachal Pradesh	27.150	92.464
1X51	Chinese Academy of Sciences	Liocichla omeiensis	China	Emei Mountain, Sichuan Province	29.52	103.33
KU15195	University of Kansas Biodiversity Institute	Liocichla phoenicea	Myanmar	Jed Lwe	26.181	98.334
KU10078	University of Kansas Biodiversity Institute	Liocichla ripponi	China	Guangxi, Diding Headwater Nature Preserve	23.122	105.963
LSUB20702	Louisiana State University Museum of Natural History	Liocichla ripponi	Captivea
LSUB37270	Louisiana State University Museum of Natural History	Liocichla ripponi	Captivea
TESRI5635	Taiwan Endemic Species Research Institute	Liocichla steerii	Taiwan	Hwalien, Hsiow-Lin, Pilu	24.182	121.308
TESRI5636	Taiwan Endemic Species Research Institute	Liocichla steerii	Taiwan	Hwalien, Hsiow-Lin, Kuanyun	24.186	121.340
TESRI5708	Taiwan Endemic Species Research Institute	Liocichla steerii	Taiwan	Nantou, Hsini, Wan Hsiang	23.619	120.928
TESRI5714	Taiwan Endemic Species Research Institute	Liocichla steerii	Captivea
TESRI5768	Taiwan Endemic Species Research Institute	Liocichla steerii	Captivea
KU15203	University of Kansas Biodiversity Institute	Leiothrix argentauris	Myanmar	Jed Lwe	26.181	98.334
KU6724	University of Kansas Biodiversity Institute	Leiothrix lutea	China	Hunan, Dongan, Dongan Shun Huang Shan National Park	26.409	111.035
KU11346	University of Kansas Biodiversity Institute	Minla ignotincta	China	Guizhou, Daozhen, Dashahe Nature Preserve	29.167	107.575
aSpecimen of captive origin with taxon determination made by examination voucher specimen plumage.

 | Show Table

DownLoad: CSV

Whole genomic DNA was extracted using a standard phenol-chloroform protocol followed by ethanol precipitation. The mitochondrial encoded cytochrome b and NADH dehydrogenase subunit 2 (ND2) were amplified for all ingroup and outgroup samples across all taxa and the mitochondrial encoded NADH dehydrogenase subunit 3 (ND3) gene was amplified from all ingroup and outgroup samples except for L. bugunorum. Conserved autosomal nuclear intron loci including introns from a magmas-like protein (12630), phosphotyrosyl phosphatase (14572), and myeloid leukemia factor 2 (23361) described in Backström et al. [34], intron 5 from transforming growth factor β2 (TGFB2, [35]) and intron 5 from β-fibrinogen (Fib5, [36]) were also amplified. In cases where primers derived from the available literature failed to amplify all samples we constructed Liocichla-specific internal primers based on sequence alignments from Liocichla and outgroup sequences.

Common sequence tags, either CS1 (5′-ACACTGACGACATGGTTCTACA) in the forward direction or CS2 (5′-TACGGTAGCAGAGACTTGGTCT) in the reverse direction, were added to all the oligonucleotide sequences in this study. Oligonucleotide tags are commonly used to facilitate the addition of adapters and barcodes to amplicons for next-generation sequencing [37]. In our study the addition of common 5′ tags allowed for CS1 and CS2 to be utilized as sequencing primers for all amplicons across all loci. A complete list of the oligonucleotide sequences used as PCR primers in this study and their properties may be found in Table 2.

Table  2.  Oligonucleotide sequences (5′-3′) used as PCR primers in this study

Oligonucleotide name	5′-CS tag + target sequence (5′ → 3′)	Locus description	Location	References
CS1-12630F	ACACTGACGACATGGTTCTACACAGCAGATCCTCAACGTCTC	Magmas-like protein	Autosomal ch: 14	[34]
CS2-12630R	TACGGTAGCAGAGACTTGGTCTCTGCAGGTAGAAGGAGCCTC	Magmas-like protein	Autosomal ch: 14	[34]
CS1-Liocichla_CS12630_F1	ACACTGACGACATGGTTCTACACAAGTGGTCGTAGTTCTGCC	Magmas-like protein	Autosomal ch: 14	This study
CS2-Liocichla_CS12630_R1	TACGGTAGCAGAGACTTGGTCTGAGATCCAGAAGGTAGAGGC	Magmas-like protein	Autosomal ch: 14	This study
CS1-14572F	ACACTGACGACATGGTTCTACACAGTAAAGAAACAGAAGTCC	Phosphotyrosyl phosphatase	Autosomal ch: 1	[34]
CS1-14572R	TACGGTAGCAGAGACTTGGTCTACTGCTGTGTGTTAGACTG	Phosphotyrosyl phosphatase	Autosomal ch: 1	[34]
CS1-Liocichla_CS14572_F1	ACACTGACGACATGGTTCTACAGTCCTCGAGACTCACATTCA	Phosphotyrosyl phosphatase	Autosomal ch: 1	This study
CS2-Liocichla_CS14572_R1	TACGGTAGCAGAGACTTGGTCTGCCATTCCTTCATAAGCTGC	Phosphotyrosyl phosphatase	Autosomal ch: 1	This study
CS1-23361F	ACACTGACGACATGGTTCTACAAAAGCTTATCAGGAGACCTC	Myeloid leukemia factor 2	Autosomal ch: 1	[34]
CS2-23361R	TACGGTAGCAGAGACTTGGTCTTTGATGTAGTCCTGCCTCTC	Myeloid leukemia factor 2	Autosomal ch: 1	[34]
CS1-Liocichla_CS23361_F1	ACACTGACGACATGGTTCTACAAGCCGCTGTCCGAGTCCCTT	Myeloid leukemia factor 2	Autosomal ch: 1	This study
CS2-Liocichla_CS23361_R1	TACGGTAGCAGAGACTTGGTCTCCCGCAGCACTTTGGCTTTGC	Myeloid leukemia factor 2	Autosomal ch: 1	This study
CS1-TGFB2F5	ACACTGACGACATGGTTCTACAGAAGCGTGCTCTAGATGCTG	Transforming growth factor β2	Autosomal ch: 3	[35]
CS2-TGFB2R6	TACGGTAGCAGAGACTTGGTCTAGGCAGCAATTATCCTGCAC	Transforming growth factor β2	Autosomal ch: 3	[35]
CS1-Liocichla_TGF_F1	ACACTGACGACATGGTTCTACATGCACACCCTCATTGTCAGACCCA	Transforming growth factor β2	Autosomal ch: 3	This study
CS2-Liocichla_TGF_R1	TACGGTAGCAGAGACTTGGTCTACAGGCAGGCAAGTCTGAGTCAC	Transforming growth factor β2	Autosomal ch: 3	This study
CS1-FIB5F	ACACTGACGACATGGTTCTACACGCCATACAGAGTATACTGTGACAT	β-fibrinogen intron 5	Autosomal ch: 4	[36]
CS2-FIB5R	TACGGTAGCAGAGACTTGGTCTGCCATCCTGGCGATTCTGAA	β-fibrinogen intron 5	Autosomal ch: 4	[36]
CS1-L5216	ACACTGACGACATGGTTCTACAGGCCCATACCCCGRAAATG	Mitochondrially encoded NADH dehydrogenase 2	Mitochondrial	[38]
CS2-H6313	TACGGTAGCAGAGACTTGGTCTACTCCTRTTTAAGGCTTTGAAGGC	Mitochondrially encoded NADH dehydrogenase 2	Mitochondrial	[38]
CS1-L10755	ACACTGACGACATGGTTCTACAGACTTCCAATCTTTAAAATCTGG	Mitochondrially encoded NADH dehydrogenase 3	Mitochondrial	[39]
CS2-H11151	TACGGTAGCAGAGACTTGGTCTGATTTGTTGAGCCGAAATCAAC	Mitochondrially encoded NADH dehydrogenase 3	Mitochondrial	[39]
CS1-L14851	ACACTGACGACATGGTTCTACACCTACTTAGGATCATTCGCCCT	Mitochondrially encoded cytochrome b	Mitochondrial	[40]
CS2-Hb745	TACGGTAGCAGAGACTTGGTCTTTTCTGGGTCTCCTAGTAGGTT	Mitochondrially encoded cytochrome b	Mitochondrial	[41]
CS1-Liocichla_cytb_F1	ACACTGACGACATGGTTCTACACATATGCCGAAACGTCCA	Mitochondrially encoded cytochrome b	Mitochondrial	This study
CS1-Liocichla_cytb_F2	ACACTGACGACATGGTTCTACACTTTCACATCGGCCGAGG	Mitochondrially encoded cytochrome b	Mitochondrial	This study
CS2-Liocichla_cytb_R1	TACGGTAGCAGAGACTTGGTCTCCTCAGAATGATATTTG	Mitochondrially encoded cytochrome b	Mitochondrial	This study
CS2-Liocichla_cytb_R2	TACGGTAGCAGAGACTTGGTCTGTCATTCTACTAGGG	Mitochondrially encoded cytochrome b	Mitochondrial	This study
CS2-Liocichla_cytb_R3	TACGGTAGCAGAGACTTGGTCTTAGTGGGTTGTTTGATCC	Mitochondrially encoded cytochrome b	Mitochondrial	This study
CS2-Liocichla_cytb_R4	TACGGTAGCAGAGACTTGGTCTGGTGTAGTAGGGGTGGAA	Mitochondrially encoded cytochrome b	Mitochondrial	This study
All oligonucleotides were synthesized by Integrated DNA Technologies (IDT) at the 25 nmole scale and purified by a standard desalting protocol. All oligonucleotides were tagged at the 5′ end with common sequence tags, either CS1 (5′-ACACTGACGACATGGTTCTACA) in the forward direction or CS2 (5′-TACGGTAGCAGAGACTTGGTCT) in the reverse direction. Locations for each locus are based on homology with the chicken (Gallus gallus) genome determined from an NCBI 'gene' search of the locus description in the case of transforming growth factor β2 and β-fibrinogen intron 5. For the remainder of the nuclear autosomal loci location was determined from the literature [34]. Internal oligonucleotides designed specifically to match conserved sequences for our samples (Liocichla and outgroups) are prefaced with 'Liocichla' in the oligonucleotide name.

 | Show Table

DownLoad: CSV

Polymerase chain reaction (PCR) was performed in 10 μL reactions with Fidelitaq™ Master Mix (Affymetrix/USB) and forward and reverse primers (0.5 μmol L^-1 each; Integrated DNA Technology) using a thermal profile of 94℃ for 4 min followed by 30 cycles of 1 min at 94℃, 1 min at 50℃, and 2 min at 72℃, and then a final extension cycle for 10 min at 72℃. The same thermal profile was utilized for all PCRs across all loci. Primers and excess dNTPs were inactivated with ExoSAP-IT^® (Affymetrix/USB) following the manufacturer's instructions. Sequencing reactions were performed on ExoSAP-IT^® treated PCR amplicons using BigDye^® v3.1 Cycle Sequencing Kit (Life Technologies). Unincorporated BigDye^® terminators were removed from sequencing reactions using the BigDye Xterminator^® kit (Life Technologies) and cleaned sequencing products were separated on a 3130 Genetic Analyzer (Life Technologies). PCR amplicons were sequenced in both directions with primers based on the common sequence tags incorporated into all amplicons (CS1 and CS2). Complementary strands were aligned and edited using Geneious v8.1.4 (Biomatters, available from http://www.geneious.com [42]). Chromatograms were inspected individually and every variation was checked for authenticity. The phase of nuclear alleles was determined computationally using the PHASE v2.1 [43] algorithms implemented in DnaSP v5 [44]. All sequences have been deposited in NCBI GenBank (Additional file 1: Table S1).

Phylogenetic analysis

We looked for evidence of recombination using the Φ_w-statistic [45] implemented in the program SplitsTree v4.10 [46]. The Φ_w-statistic distinguishes recombination from recurrent mutation. We tested for selection using the Hudson-Kreitman-Aguade (HKA) test [47] implemented in DnaSP v5 [44]. We selected the best-fit model of nucleotide substitution for each locus using Akaike information criteria (AIC) in the program jModeltest2 v2.1.7 [48]. We used PhyML v3.0 [49, 50] implemented in Geneious v8.1.4 [42] to generate maximum likelihood gene trees and tested the molecular clock hypothesis for each locus using a likelihood ratio test in PAUP^* v4.0b10 [51]. Maximum likelihood gene tree analyses were performed using a concatenated dataset of all three mtDNA genes as well as for each nuclear locus. Bootstrap support was determined with 1, 000 replicates.

We combined all sequence information into a locus partitioned ^*BEAST mutlilocus species tree analysis implemented in the BEAST v2.3.0 software platform [52, 53]. ^*BEAST uses a Bayesian Markov chain Monte Carlo (MCMC) method to jointly estimate multiple gene trees embedded within a shared species tree. Individual sequences were grouped by taxon and each taxon was used as a separate operational taxonomic unit (OTU) in all species tree analyses; L. bugunorum and Minla ignotincta were excluded because each lacked sequence data for at least one locus. Because of their uniparental inheritance and lack of recombination, mitochondrial genes were constrained to the same tree model but given gene-specific substitution and clock models. Substitution models were determined for each locus using the bModelTest package [54] implemented in BEAST v2.3.0. The molecular clock hypothesis was rejected for all three mtDNA genes as well as Fib5, so these loci were assigned uncorrelated lognormal relaxed clock models [55]; strict clock models were used for all other loci. For species tree priors, we used a Yule tree prior and a piecewise linear and constant root population size model.

Divergence dates were estimated using a molecular clock approach. Lineage substitution rates were estimated for all loci in BEAST v2.3.0 based on lognormal prior distributions. The mitochondrial cytochrome b substitution rate was estimated based on a lognormal prior distribution with a mean and initial value of 0.0105 substitutions/site/lineage/Ma [56] and a standard deviation of 0.15 (95% interquartile range of 0.008-0.0139). This distribution is overlapping with other estimates for cytochrome b substitution rates in birds [56-58]. All other loci were estimated based on more broadly distributed lognormal prior distributions for substitution rate. ND2 and ND3 substitution rates were estimated based on a lognormal distribution with a mean of 0.0105 substitutions/site/lineage/Ma and a standard deviation of 0.75 (95% interquartile range of 0.002-0.0345). All nuclear loci substitution rates were estimated based on a lognormal prior distribution with a mean of one-tenth of the mean substitution rate of the mitochondrial loci (0.00105 substitutions/site/lineage/Ma) and a standard deviation of 0.75 (95% interquartile range of 0.0002-0.0035).

MCMCs were run for 100 million generations (sampling every 10, 000 generations and discarding the first 10% as burn-in). Four independent MCMCs were run under these conditions and log files and species tree files from each run were combined into a single log file and a single species tree file using the logCombiner software included with BEAST v2.3.0. Convergence was assessed in the program Tracer v1.6.0 [59].

We also conducted two additional analyses using ^*BEAST in BEAST v2.3.0. First, we assessed the influence of codon partitioning on both species tree topology and species tree height by partitioning each of the three mitochondrial coding genes by all three codon positions and running the same analysis in four replicate runs. Third position codon partitions were given lognormal prior distributions for lineage substitution rate with the third position for cytochrome b with a mean and initial value of 0.0105 substitutions/site/lineage/Ma and a standard deviation of 0.15 and ND2 and ND3 third position sites with a mean value of 0.0105 substitutions/site/lineage/Ma and standard deviation of 0.75. First and second codon positions were given lognormal prior distributions with a mean of one-tenth the third position substitution rate and a standard deviation of 0.75. All other parameters were the same as those analyses with partitioning only by locus. We evaluated models that used both partitioning by locus and codon against those with data partitioned only by locus in relation to species tree height by comparing marginal likelihoods using the Akaike's information criterion (AIC) from the model comparison program in Tracer v1.6.0 [59, 60]. Additionally, to confirm that mtDNA data were not driving the inferred species tree topology, we conducted a separate ^*BEAST analysis in BEAST v2.3.0 under identical conditions, save for the exclusion of the three mtDNA genes. All relevant xml input files and output files are available on Dryad (http://dx.doi.org/10.5061/dryad.413f4).

Results

We obtained a total of 1, 733 bp of mtDNA sequence data and 1, 921 bp of nuclear intron sequence data across 11 ingroup samples representing all five Liocichla species (L. omeiensis, L. bugunorum, L. steerii, L. phoenicea and L. ripponi) and three outgroup samples from three outgroup species (Leiothrix argentaurus, Leiothrix lutea and Minla ignotincta, Table 3). Data from all loci were obtained from at least one individual sample across all taxa except L. bugunorum, for which we were only able to obtain sequences from the cytochrome b (KJ456321.1) and ND2 (KJ455478.1) genes [33], and M. ignotincta for which we were missing data for one nuclear intron loci requiring us to exclude this taxon from the ^*BEAST species tree analysis. The 14572 intron had a 4-bp insertion/deletion (indel) within the ingroup that was identified by eye and removed prior to analyses. Four additional indels were identified in the outgroup taxa and also removed prior to analyses. 42 base pairs of tRNA sequence were excluded from ND3 for the ^*BEAST analysis to facilitate codon partitioning. No locus exhibited evidence of recombination.

Table  3.  Descriptive statistics for the eight molecular loci used to construct phylogenetic relationships within the genus Liocichla

Locus	Inheritance	# Sequences	Length (bp)	Alleles	s	PI sites	π (×10-3)	Model
cytb	Mitochondrial	11	400	10	68	51	72.95	TPM2uf + I
ND2	Mitochondrial	11	949	11	204	133	84.99	TrN + G
ND3	Mitochondrial	10	384	5	64	42	69.97	TPM3uf + G
FIB5	Autosomal	20	511	6	11	8	7.82	TIM3 + G
TGFB2	Autosomal	20	359	9	22	14	17.83	TrNef + I
12630	Autosomal	16	203	10	13	8	19.95	TPM2uf
14572	Autosomal	20	451	10	14	5	6.44	HKY + G
23361	Autosomal	16	397	5	11	6	8.35	HKY + I
The summary includes the assumed inheritance, the number of sequences analyzed, the length of the sequences (in base pairs), the number of alleles, the number of segregating sites (s), the number of parsimony-informative sites (PI sites), nucleotide diversity (π), and the estimated best-fit nucleotide substitution model. Statistics are drawn from ingroup sequences only, whereas best-fit nucleotide substitution models were estimated with outgroups.

 | Show Table

DownLoad: CSV

Maximum likelihood analysis for mtDNA sequences and MCMC Bayesian phylogenetic analyses in ^*BEAST for both a combined nuDNA and mtDNA dataset and a nuDNA only dataset recovered the same strongly supported topology, thereby providing robust support for a single phylogenetic hypothesis for the genus Liocichla. Individual gene trees derived from analysis of the nuclear introns showed low resolution (Fig. 2). Comparatively poor resolution in nuDNA gene trees is expected given the larger effective population sizes of nuclear genes versus mitochondrial genes [61]. Maximum-likelihood analysis of the concatenated mtDNA dataset resulted in a well-resolved topology that placed L. steerii as the sister taxon to all remaining extant Liocichla. L. omeiensis was recovered as sister to L. bugunorum and L. phoenicea was recovered as sister to L. ripponi (Fig. 3).

Figure  2.  nuDNA maximum likelihood gene trees. Nuclear intron maximum likelihood gene trees generated with PhyML v3.0 for nuDNA sequence data for 12630 (a), 14572 (b), 23361 (c), Fib5 (d) and TGFB2 (e). Each sequence includes the taxon, voucher institution and number and is followed by a or b to denote the phased alleles derived from a single diploid sequence. Branch support with bootstrap values greater than 90 are shown.

DownLoad: Full-Size Img PowerPoint

Figure  3.  mtDNA maximum likelihood gene tree. Maximum likelihood mitochondrial gene tree based on concatenated sequence data from cytochrome b, ND2 and ND3 (outgroup taxa not shown).

DownLoad: Full-Size Img PowerPoint

The ^*BEAST species tree analysis resulted in a strongly supported topology consistent with the mtDNA maximum likelihood gene tree topology (Fig. 3) with posterior branch support of 0.98-1.0 (Fig. 4). The same topology was produced with only moderately lower posterior branch support (0.85-1.0) when mtDNA data were excluded from the analysis (Fig. 5). A ^*BEAST species tree analysis in BEAST v2.3.0 with data partitioned by both locus and codon position resulted in the same topology with similar levels of posterior branch support (0.99-1.0) but with comparatively older divergence times relative to both an analysis of our data with partitioning only by locus and divergence estimates of Liocichla from prior studies [11]. However, a model that is partitioned by both locus and codon was shown to be a significantly worse fit with regards to species tree height (age) than a model with partitioning by locus alone.

Figure  4.  Combined mitochondrial and nuclear DNA species tree. Species tree analysis using the program ^*BEAST in BEAST v2.3.0 and combined mtDNA and nuDNA sequence data partitioned by locus. Branches labeled with posterior branch support values.

DownLoad: Full-Size Img PowerPoint

Figure  5.  Nuclear DNA species tree. Species tree analysis using the program ^*BEAST in BEAST v2.3.0 using only nuDNA sequence data partitioned by locus. Branches labeled with posterior branch support values.

DownLoad: Full-Size Img PowerPoint

Divergence time estimation based on a ^*BEAST analysis incorporating partitioning by locus only indicated that L. steerii began its divergence from the rest of the extant Liocichla approximately 8.94 Ma [5.55-12.87 Ma highest posterior density (HPD)]. The next split between the L. omeiensis and L. phoenicea/L. ripponi clade occurred approximately 6.06 Ma (3.06-9.18 Ma HPD). The split between L. phoenicea and L. ripponi occurred approximately 0.88 Ma (0.07-1.88 Ma HPD).

Estimated samples sizes (ESS) from the combined ^*BEAST locus partitioned runs for all posterior mean lineage substitution rates for relaxed lognormal clock models and posterior lineage substitution rates for strict clock models were over 10, 000 for all loci except for ND2 (ESS_ucldmean = 4695). ESS values from the combined ^*BEAST runs for the species coalescent, species tree height and Yule species tree model were all > 6, 000 with consistent traces for each.

Discussion

Phylogenetic relationships

Prior work supports the existence of a monophyletic Liocichla closely allied with the genera Actinodura, Heterophasia, Leiothrix, Minla but none of these prior studies had sufficient taxon sampling to investigate the relationships within Liocichla [6, 10, 11]. The origin of the Liocichla separates an island endemic (L. steerii) from a continental lineage containing all other taxa (Figs. 3, 4, 5). L. bugunorom, L. omeiensis and L. steerii share several plumage characters and Collar and Robson [1] have hypothesized that these taxa form a "superspecies", however, geography in this case seems to be a better indicator of evolutionary history than morphology. The plumage similarity between L. bugunorum and L. steerii is particularly striking, and several features, including the yellow lores, yellow-tipped undertail coverts, and overall green coloration, are uniquely shared by these two Liocichla species [1, 25, 62]. Our phylogenetic analysis, however, rejects a close relationship between L. steerii and either L. bugunorom or L. omeiensis and suggests the observed plumage similarities are homoplasious, possibly as a result of retained ancestral characters. Similar results have been described among other Asian birds [63-65] and discordance between morphology and molecular based phylogenies is widespread within the Sylvioidea [4]. These studies highlight the labile nature of avian plumage characters in general [66-68].

Biogeography

The climatic and geologic history of Asia has given rise to a complex and dynamic landscape with a myriad of potential isolating barriers and ecological shifts that wax and wane over time. The collision of the Indian subcontinent with Asia and the resulting rise of the Himalayas, episodic climate change creating shifting precipitation patterns and periodic land bridge connections resulting from rising and falling sea levels [69-72] have all contributed to the biogeographical patterns found among Asian species [12, 31, 33, 73-78]. These events also likely played critical roles in the diversification of Liocichla in Asia.

The island of Taiwan in particular provides an important testing ground for hypotheses on the assemblage of biotas on continental islands. The phylogenetic position of L. steerii hints at a comparatively early divergence on the island of Taiwan. Alternatively, the early divergence of L. steerii may be illusory due to our inability to include data from one or more extinct sister lineages. Evidence for a sister relationship between a Taiwan endemic and an Indochinese endemic in Vietnam has been found for the Green-backed Tit (Parus monticolus). Wang et al. [79] postulated that these taxa represent the relic populations of a wider radiation that connected populations on Taiwan with those in Mainland Asia. This hypothesis would predict that for those Taiwan endemics for which there are extant lineages in the region from Southeast China to Northern Indochina that these lineages should be more closely related to Taiwan endemics than to other continental lineages. This is indeed the case for other widespread babbler radiations in Asia with an endemic taxon on the island of Taiwan. In Alcippe morrisonia the endemic taxon on Taiwan is sister to a clade containing taxa from Hainan and Fujian [80, 81] and in the widespread Pomatorhinus ruficollis complex the endemic P. musicus on Taiwan forms the sister lineage to a lineage consisting of P. nigrostellatus in South China, Hainan and Northern Indochina and P. stridulus in Eastern China [63, 65].

Molecular clock-based estimates of divergence times, while admittedly fraught with assumptions [57, 82-84], provide useful benchmarks in broad comparisons both among species and with the dating of key geological and climatic events. The divergence of L. steerii of 8.94 Ma (5.55-12.87 Ma HPD) may precede the Pliocene/Miocene boundary and is either earlier or overlapping with only the oldest estimates of the timing of Taiwan's uplift approximately 6 Ma [85-87]. The divergence of L. steerii is also either older or overlapping only with the older end of the distribution for divergences estimated for other montane forest taxa in Taiwan including birds (3.51 Ma, Alcippe morrisonia, [80]; 3.25 Ma, Regulus goodfellowi, [88]; 1.7 Ma, Carpodacus formosanus, [89]; 2.6 Ma, Carpodacus formosanus, [90]; 0.18 Ma, Parus monticolus, [79]; 3.1 Ma, Sittiparus varius, [91]), trees (3.23-3.41 Ma, Taiwania cryptomerioides, [92]), amphibians (2.22 Ma, Sylvirana latouchii, [93]; 1.25 Ma, Limnonectes fujianensis, [94]) and insects (0.3 Ma, Euphaea formosa, [95]). L. steerii's comparatively early appearance on Taiwan may in part coincide with the early appearance of subtropical forests suitable for subtropical Oriental forest flora and fauna. Some of the later diverging endemic taxa on Taiwan are higher elevation species with a Palearctic origin whose settlement on Taiwan was likely postponed until the island had undergone sufficient uplift to foster more temperate biomes.

These data suggest that L. steerii may either be among the oldest endemic avian species in Taiwan, potentially isolated by ecological barriers that predate some currently recognized estimates for the uplift of the island, or like P. musicus [63, 65] and A. morrisonia [80, 81], may have diverged later from ancestral populations in Eastern and Southern China who in the case of L. steerii have gone extinct. The comparatively patchy distribution of L. ripponi in Southern and Eastern China suggests that the environment in this region may have been less than ideal for Liocichla.

The extant continental Liocichla consist of four distinct evolutionary lineages (L. phoenicea, L. ripponi, L. bugunorum, L. omeiensis). The present study supports a L. phoenicea/L. ripponi lineage sister to a L. bugunorum/L. omeiensis lineage (Figs. 3, 4). We were not able to incorporate sequences from L. bugunorum into a multilocus species tree, however, mitochondrial DNA supports a sister relationship between L. bugunorum and L. omeiensis with the Hengduan Mountains as a potential isolating barrier between these two species. Diversification of Liocichla species in mainland Asia follows many of the patterns observed in other avian taxa in Asia. Southcentral China and Southeast Asia are in particular pivotal centers of origin for biodiversity of Asia, particularly for passerine birds [31, 32, 74], and these regions encompass the probable ancestral Miocene distribution for the Liocichla [10, 11, 31]. Only L. bugunorum and L. phoenicea exhibit Himalayan distributions and each has a sister taxon with a different non-Himalayan distribution (Central China and South China/Northern Indochina respectively). These data suggest that the Himalayan Liocichla may have independently originated from dispersal followed by vicariance between dispersed Himalayan populations and ancestral populations in Central China and Southeast Asia. Dispersal into the Himalayas followed by vicariant divergence has also been supported by phylogenetic and biogeographic analyses for other passerine birds in the region [31, 74].

Geological and climatic change likely created new opportunities for range expansion, niche filling and subsequent diversification of the Liocichla into emerging montane forest environments the Eastern Himalaya, Central China and East Asia [31, 33, 96, 97]. The L. phoenicea/L. ripponi lineage identified in the present study provides fodder for future work. Together with diagnostic plumage differences [23], our results support species status for these taxa, recognized by previous authors as subspecies [29]. A Pleistocene divergence time for L. phoenicea and L. ripponi roughly coincides with those for some montane forest birds in the region [80] but is later than others [74]. The Mekong-Salween Divide has been proposed as a biogeographic barrier for the region [98, 99]. The divide between the parapatric distributions of L. phoenicea and L. ripponi appears to be west of the traditional Mekong-Salween Divide approximately at the Myitnge River and Myanmar's Eastern highlands [1] but additional data on the distribution of these taxa are needed. Also, the fragmented range of both L. phoenicea and L. ripponi across South and Southeast Asia may be indicative of further phylogeographic structure in these lineages that may be uncovered with additional sampling.

Conclusions

This study when taken in the context of other phylogenetic and biogeographic work in the region poses a plausible hypothesis for the evolution of Liocichla in Asia. The Liocichla radiation likely was initiated sometime between the mid-Miocene to the Miocene-Pliocene boundary. Changing climate during this period likely opened new ecological niches to Southeast Asian immigrants resulting in speciation events across a broad region from East Asia to South and Central China to the Eastern Himalayas [72, 74, 76, 78, 97].

An eastern expansion of Liocichla ended at the island of Taiwan where L. steerii ultimately became isolated, but this expansion left no extant lineages of Liocichla in Eastern and Southeastern China. Alcippe morrisonia also has an endemic taxon on the island of Taiwan but unlike L. steerii, A. morrisonia on Taiwan has extant sister lineages in Eastern and Southeastern China. However, the extant populations of A. morrisonia in Eastern and Southeastern China likely underwent significant demographic fluctuations compared to those in Central and Western China [80]. Periodic climatic change occurring from the Miocene-Pliocene boundary to the Pleistocene may have had comparatively greater effects on the ecosystems of Eastern and Southeastern China causing population bottlenecks in some taxa (A. morrisonia) [80] and extinction in others (Liocichla).

The history on the western side of the Liocichla radiation is more complex and likely involved two independent instances of range expansion into the Eastern Himalayas followed by speciation. This pattern is consistent with the conclusions drawn from studies of other Asian birds with taxa in the Eastern Himalayas [31, 74]. The most recent split in the Liocichla appears to have occurred during the Pleistocene between the parapatrically distributed L. phoenicea and L. ripponi although it is unclear as to whether or not this represents an example of parapatric speciation or secondary contact.

Future studies involving ecological niche modeling for the genus, reconstruction of ancestral distributions and more broadly sampled molecular-based studies of phylogeography and historical demography for the more widespread but fragmented L. phoenicea and L. ripponi would be valuable contributions in further testing of the hypotheses generated from this study. These studies however would all require additional sampling and field data especially for those taxa with lingering questions about range limits.

Availability of supporting data

The data set supporting the results of this article is available in the Dryad repository, http://dx.doi.org/10.5061/dryad.413f4.

Additional file

Additional file 1: Table S1. Genbank accession numbers for all ingroup and outgroup sequences used in this study

Authors' contributions

HLM conceived the study, participated in generating the molecular data in the lab, designed novel primers, conducted sequence validation and alignment, conducted phylogenetic analyses and led efforts to draft the manuscript. BDM conducted phylogenetic analyses and helped draft the manuscript. DTT provided the sequence data for Liocichla bugunorum and helped draft the manuscript. CTY acquired samples for Liocichla steerii, participated in the coordination of the study and helped draft the manuscript. LNM and KNM participated in generating the molecular data in the lab and helped with sequence validation and alignment. FL acquired a sample for Liocichla omeiensis, participated in coordination of the study and helped draft the manuscript. All authors read and approved the final manuscript.

Acknowledgements

Support was provided by a National Science Foundation (NSF) Major Research Instrumentation award (DBI-0821703) to HLM, Roger Fry and Len Weakley through the Helen B. Vogel Charitable Trust and the Mary Jane Helms Charitable Trust to HLM and Deutsche Forschungsgemeinschaft to DTT (Ti 679/1-1). Thanks to Jamie Allison and the biotechnology program at Loveland High School for identifying students to participate in the lab component of this research. Thanks to Nadja Spitzer for translating one of our references from German to English. Thanks to Mark Robbins and Robert Moyle of the Kansas University Biodiversity Institute and Donna Dittman and Fredrick Sheldon of the Louisiana State University Museum of Natural Science for grants of tissue samples. Thanks to Per Alström, Shou-Hsien Li and two anonymous reviewers for valuable comments on this manuscript.

Compliance with ethical guidelines

Competing interests

The authors declare that they have no competing interests.