Accepted Manuscript Short Communication A reinvestigation of phylogeny and divergence times of Hynobiidae (Amphibia, Caudata) based on 29 nuclear genes Meng-Yun Chen, Rong-Li Mao, Dan Liang, Masaki Kuro-o, Xiao-Mao Zeng, Peng Zhang PII: DOI: Reference:

S1055-7903(14)00361-3 http://dx.doi.org/10.1016/j.ympev.2014.10.010 YMPEV 5050

To appear in:

Molecular Phylogenetics and Evolution

Received Date: Revised Date: Accepted Date:

21 August 2014 13 October 2014 14 October 2014

Please cite this article as: Chen, M-Y., Mao, R-L., Liang, D., Kuro-o, M., Zeng, X-M., Zhang, P., A reinvestigation of phylogeny and divergence times of Hynobiidae (Amphibia, Caudata) based on 29 nuclear genes, Molecular Phylogenetics and Evolution (2014), doi: http://dx.doi.org/10.1016/j.ympev.2014.10.010

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.

Title: A reinvestigation of phylogeny and divergence times of Hynobiidae (Amphibia, Caudata) based on 29 nuclear genes

Authors: Meng-Yun Chen1#, Rong-Li Mao1#, Dan Liang1, Masaki Kuro-o2, Xiao-Mao Zeng3, Peng Zhang1*

Address: 1

State Key Laboratory of Biocontrol, School of Life Sciences, Sun Yat-Sen

University, Guangzhou 510006, China 2

Department of Biology, Faculty of Agriculture and Life Science, Hirosaki

University, Hirosaki, Japan 3

Chengdu Institute of Biology, Chinese Academy of Sciences, Chengdu 610041,

China

*

Corresponding author:

Peng Zhang. #434, School of Life Sciences, Sun Yat-Sen University, Higher Education Mega Center, Guangzhou 510006, China; Email: [email protected]

# these authors equally contribute to this work.

Abstract Although several recent studies have investigated the major phylogenetic relationships within Hynobiidae, their evolutionary history remains partially resolved and the phylogenetic positions of some genera, particularly Pachyhynobius and Salamandrella are still disputed. Notably, previous studies relied primarily on mitochondrial DNA data and concatenated analyses; thus, a new investigation based on multiple nuclear genes and species-tree inference is needed. Here, we provide an in-depth phylogenetic analysis, based on 29 nuclear genes comprising 29,232 bp of data from a comprehensive taxonomic sampling (24 hynobiids and 7 outgroups), using both concatenated and species-tree methods. Our results robustly resolved most genus-level relationships within Hynobiidae, including the placement of Salamandrella as the sister group to a clade containing Batrachuperus, Liua and Pseudohynobius, and the placement of Pachyhynobius as the sister group to a clade containing all hynobiids excluding Onychodactylus, Paradactylodon and Ranodon. Time estimates based on our data suggest that the major group of living hynobiids (excluding Onychodactylus) originated approximately 40 Ma, ~50% younger than estimates from mtDNA data (62.5 Ma) but 10% older than estimates from three nuclear genes (36 Ma). Our results highlight the benefits of using a large number of nuclear loci to infer both phylogeny and time for relatively old lineages.

Keywords: molecular dating; salamander, Pachyhynobius, Salamandrella

1. Introduction The family Hynobiidae is the third largest taxonomic family of living salamanders, including 9 genera and 59 currently recognized species (AmphibiaWeb 2014). Morphologically, Hynobiidae, together with the family Cryptobranchidae, possess some primitive traits such as external fertilization, an angular bone in the lower jaw, and large numbers of microchromosomes (Regal 1966; Edwards 1976). Recent molecular analyses also separate Hynobiidae from other extant salamanders early in salamander phylogeny (Roelants et al. 2007; Pyron & Wiens 2011; Shen et al. 2013). Five different studies have investigated phylogenetic relationships on a broad scale within the Hynobiidae. Using complete mitochondrial genome data from 16 hynobiid species, Zhang et al. (2006) presented a well resolved phylogeny for most relationships within the family. However, another analysis also based on mitochondrial genomes with an improved taxon sampling modified the position of Pachyhynobius, although not conclusively (Peng et al. 2010). More recently, Zheng et al. (2011) used mtDNA genome data and three nuclear genes to investigate phylogeny and divergence times for Hynobiidae, including 25 hynobiid species with representatives of all genera. In that study, the phylogenetic tree inferred from mitochondrial genomes was nearly identical to that of Zhang et al. (2006), whereas different placements of Pachyhynobius and Salamandrella were obtained using nuclear genes. As part of a phylogenetic study of all amphibians, Pyron and Wiens (2011) also estimated phylogenetic relationships for 49 hynobiid species using DNA

sequence data (primarily mtDNA) collected from GenBank. This study yielded some different intergeneric relationships for Hynobiidae, but generally with low bootstrap support. Finally, using two mitochondrial genic regions but dense taxon sampling, Weisrock et al. (2013) found many discordant branches with low support among their different analyses and considered the major intergeneric relationships of Hynobiidae a polytomy. Current estimates of the deeper regions of hynobiid phylogenetic history rely primarily on mitochondrial DNA data and concatenated analyses. It is well known that mtDNA as a whole is a single unit of genetic transmission; thus, mtDNA-based studies would better be cross validated by multiple independent nuclear genes. Furthermore, It is now realized that traditional data concatenation, which combines alignments in a supermatrix, sometimes can produce incorrect results with high confidence (Edwards 2009; Lemmon & Lemmon 2013). Therefore, analyses based on multiple nuclear genes and non-concatenated framework are necessary to resolve hynobiid phylogeny (Weisrock et al. 2013). In addition, previous dating analyses for hynobiid evolution using mitochondrial genomes (Zhang et al. 2006) or a few of nuclear genes (Zheng et al. 2011) have also yielded conflicting results. Times derived from the mitochondrial data are much older than those estimated from nuclear data when the reference point for calibration is a fossil-based estimate of an event preceding the phylogenetic divergences being studied. Recent empirical studies have consistently shown that slowly evolving markers such as nuclear exons surpass mitochondrial data for dating divergences

exceeding 15 Mya (Zheng et al. 2011; Near et al. 2012). However, because nuclear genes are relatively conservative and contain fewer signal sites, more loci are needed to calculate reliable divergence times. Here, we reinvestigate the phylogenetic relationships and divergence times of Hynobiidae, using both concatenated and non-concatenated analyses with 29 nuclear genes. The resulting time estimates are compared to previous analyses based on different datasets. We expect to resolve robust phylogenetic relationships among the major clades of Hynobiidae and to obtain more accurate divergence-time estimation to improve our understanding of the evolution of this family.

2. Materials and methods 2.1 Taxon sampling and DNA sequencing 24 hynobiid species were sampled in this study, representing all nine currently recognized genera of Hynobiidae (AmphidiaWeb 2014). Another seven non-hynobiid salamander species were used as outgroups, covering six salamander families other than Hynobiidae. Details for taxon sampling are presented in Table S1. Total genomic DNA was extracted from muscle or liver tissue using standard salt extraction methods. The PCR primers and protocol to amplify the 29 nuclear loci are from a nuclear marker toolkit recently described by Shen et al. (2013). This toolkit includes 102 nuclear markers (all protein-coding sequences) that can be universally applied in most vertebrate groups and thus fits our research. PCR products were purified by ExoSAP (USB) treatment and bidirectionally sequenced with Sanger sequencing.

2.2 Phylogenetic analyses Sequences of each nuclear gene were separately aligned with PRANK (Löytynoja & Goldman 2008) at default settings according to their translated amino acid sequences (-translate). Ambiguous alignment regions were trimmed by using Gblocks (Castresana 2000) with type of sequences set to codons (-t=c) and half gaps allowed (-b5=h); otherwise, default settings were assumed. The concatenated 29-locus dataset was analyzed with both maximum likelihood (ML) and Bayesian inference (BI) approaches. The best-fit partitioning scheme and corresponding nucleotide substitution models were selected using the Bayesian information criterion (BIC) implemented in PartitionFinder (Lanfear et al. 2012). The ML analysis was conducted using RAxML version 7.2.6 (Stamatakis 2006). Rapid bootstrap analysis (-f a) with 500 replicates was performed to estimate branch support. The topological test for comparing alternative hypotheses was done by using CONSEL (Shimodaira and Hasegawa 2001) with site-wise log-likelihoods calculated by RAxML. The Bayesian analysis was conducted using MrBayes 3.2 (Ronquist et al. 2012). Four Markov chains were run with one cold chain and three heated chains (temperature = 0.1) . Two independent runs were performed for 20 million generations, and trees were sampled every 500 generations. The stationarity of the likelihood scores of sampled trees was assessed using Tracer version 1.5. The first 25% of generations were discarded as burn-in. Two different approaches were used to estimate species trees under a coalescent

framework, a pseudo-likelihood method implemented in the program MP-EST (Liu et al. 2010) and a Bayesian MCMC method implemented in the program BEAST version 1.7.5 (Drummond & Rambaut 2007). Both methods constructed the species tree based on gene trees estimated from the 29 nuclear loci. In BEAST analyses, each locus was assigned a separate GTR+G+I model, and the substitution models, clock models, and tree models for the 29 nuclear loci were unlinked. A relaxed uncorrelated lognormal clock and a Yule speciation prior were assumed. Two independent MCMC runs were conducted for 500 million generations, sampling every 50,000 generations. The convergence of the independent runs was assessed in Tracer version 1.5. The first 200 million generations were discarded as burn-in. In MP-EST analyses, gene trees were reconstructed for each locus under the GTR+G+I model using PHYML 3.0 (Guindon et al. 2010). The gene trees generated from PHYML were used as input to generate MP-EST trees using the program MP-EST. The statistical support was estimated with 500 bootstrap replicates. Bayesian concordance analysis (Ané et al. 2007) was conducted to assess the level of concordance in phylogenetic analyses across different loci. We performed the analysis using the program BUCKy version 1.4.0(Ané et al. 2007). All the individual gene trees needed in BUCKy analyses were generated in MrBayes 3.2. The BUCKy program was run for 10 million generations after 1 million generations of burn-in. Two independent replicate runs were performed. We tried a range of Dirichlet priors (α= 0.1, 1, 2.5, 10) in the BUCKy analyses and found no effect of the prior on the analyses.

2.3 Estimating divergence time Divergence times were estimated using two relaxed-clock based programs: MultiDivTime (Thorne and Kishino 2002) and MCMCTREE (Dos Reis et al. 2014). In both analyses, the ingroup root is the Cryptobranchidae-Hynobiidae split, constrained at 155 ± 10 million years ago (Mya). This was based on the earliest known cryptobranchid fossil Chunerpeton tianyiensis from Late Jurassic (Gao and Shubin 2003). In the MultiDivTime analysis, optimized branch lengths with their variance–covariance matrices were estimated for each gene with an F84 + G model. The priors for the mean (rttm) and standard deviation (rttmsd) of the ingroup root age were set to 1.55 and 0.1, respectively. The prior mean (rtrate) and standard deviation (rtratesd) for the gamma distribution describing the rate at the root node were both set to 0.03669. These values were based on the median of the substitution path lengths between the root and each terminal, divided by rttm (as suggested by the authors). The autocorrelation parameter prior (brownmean) and its SD (brownsd) were set to 1.29, such that brownmean multiplied by the rttm prior (1.55) equals 2.0. After an initial burn-in period of 1,000,000 cycles, MCMC chains were run for 4,000,000 cycles, with sampling intervals of every 200 cycles. Two independent runs were performed to examine whether similar results were observed. In the MCMCTREE analysis, the approximate likelihood method was used for divergence time estimation (see MCMCTREE tutorial), with a separate GTR+G

model used for each locus. Two models of molecular clock were used: the independent rates (clock=2) and the autocorrelated rates (clock=3). The root age is set at B(145, 165), using the soft-bound strategy with 2.5% tail probabilities above or below the limits. The first 500,000 generations were discarded as burn-in, and the Markov chain was sampled every 200 generations until 20,000 samples were collected. Two runs were performed to ensure convergence of the Markov chain.

3. Results and discussion 3.1 Phylogeny of hynobiids A total of 639 new sequences were generated for this study. All obtained sequences have been deposited in GenBank under accession numbers KJ715237-KJ715875 (Table S1). The supermatrix of 29 combined nuclear genes is 98.9% complete for all 31 salamander taxa and 99.0% complete for 24 hynobiid taxa. The aligned lengths of 29 nuclear genes range from 549 to 1503 bp (mean = 1008). Gene names, alignment lengths, taxon coverages, variable-site information, and selected substitution models are shown for each locus in Table 1. The concatenated 29-gene dataset comprises 29,232 base pairs (bp) and displays no apparent substitutional saturation within all 31 included salamander species (Fig. S1). Both maximum-likelihood and Bayesian analyses of the 29 nuclear gene dataset produced identical topologies. These Bayesian and ML trees had similar branch lengths and levels of support; a total of 20 (of 23) internal nodes in the ML tree

had >70% bootstrap support, and all internal nodes in the Bayesian tree had posterior probabilities > 0.95 (Fig. 1). The two species-tree analyses (MP-EST and BEAST) produced topologies and branch support highly similar to those from the concatenated analyses (Fig. 1). Many branches in the final tree received high concordance factors (CFs) in the Bayesian concordance analysis (BUCKy), suggesting that the final tree is supported by most of the 29 nuclear genes (Fig. S2). For each genus represented by two or more species, monophyly of the genus is strongly supported (Fig. 1). Nearly all the intergeneric relationships within Hynobiidae are robustly resolved in this study. The genus Onychodactylus is resolved as the sister lineage of all other hynobiids, and Ranodon is placed in a clade with Paradactylodon, as recovered by all recent molecular studies (Zhang et al., 2006; Peng et al., 2010; Zheng et al., 2011; Pyron and Wiens 2011; Weisrock et al., 2013). The placement for the monotypic genus Pachyhynobius is one of the most compelling questions in hynobiid phylogenetics. The most likely position of Pachyhynobius inferred from mitochondrial genomes is as the sister taxon to all remaining hynobiids except Onychodactylus, Ranodon and Paradactylodon (Zhang et al., 2006); however, this result was strongly supported only in their unpartitioned Bayesian analysis. By using three nuclear genes, Zheng et al. (2011) recovered Pachyhynobius as the sister taxon to a clade excluding Onychodactylus, Ranodon, Paradactylodon and Salamandrella, but without statistical support. Pyron and Wiens (2011) found Pachyhynobius as the sister taxon to Salamandrella, whereas the branch support was low and there existed extensive missing data in their dataset. Weisrock et

al. (2013) had not resolved the position of Pachyhynobius in their analyses with two mitochondrial regions. Our concatenated analyses based on 29 nuclear genes strongly support Pachyhynobius as the sister lineage of all other hynobiids excluding Onychodactylus, Paradactylodon and Ranodon (MLBS > 95%, BPP = 1.00; Fig. 1), in agreement with the previous hypothesis based on complete mitochondrial genomes (Zhang et al., 2006). This result is also strongly supported in the two species-tree analyses, which accounts for potential discordance among individual gene histories (MP-EST BS > 95%, BEAST PP = 1.00; Fig. 1). Importantly, all other possible hypotheses about the placement of Pachyhynobius can be rejected by both KH and AU tests (Table S2). Hynobius, Salamandrella, Batrachuperus, Liua and Pseudohynobius form a well-supported clade. Within this clade, Salamandrella is the sister taxon of another well-supported clade containing Batrachuperus, Liua and Pseudohynobius (Fig. 1). This relationship is strongly supported in the ML, Bayesian and BEAST analyses (MLBS = 94%, BPP = 1.0, BEAST-PP = 1.0; Fig. 1) and moderately supported in the MP-EST analysis (BS = 64%; Fig. 1). In addition, other hypotheses related to the position of Salamandrella can be rejected by both KH and AU tests (Table S2). This result is consistent with previous molecular studies based on complete mitochondrial genomes (Zhang et al., 2006; Peng et al., 2010; Zheng et al., 2011), but contrary to the studies based on a few nuclear genes (Pyron and Wiens 2011; Zheng et al., 2011), suggesting that a substantial number of nuclear genes are needed to resolve phylogenetic hypotheses.

Among the three montane hynobiid genera Batrachuperus, Liua and Pseudohynobius, the aquatic Liua is closer to the terrestrial Pseudohynobius rather to the aquatic Batrachuperus, in congruence with most previous molecular studies (Zhang et al., 2006; Peng et al., 2010; Zheng et al., 2011; Pyron and Wiens 2011; Weisrock et al., 2013). However, none of the studies (including ours) have robustly resolved the relationships among these three genera. The internal branch separating these genera is extremely short (node G; Fig. 1), suggesting that these ancient cladistic events occurred almost simultaneously. Although not fully resolved, the degree of resolution to this part is nevertheless encouraging, and holds promise that using increasingly large datasets, especially from additional nuclear protein-coding loci, may resolve this problem.

3.2 Divergence times of hynobiid evolution In this study, we used 29 nuclear genes and three relaxed-clock dating algorithms, MultiDivTime (correlated-rate) and MCMCTREE (both independent-rate and correlated-rate models) to calculate divergence times for hynobiids (Table 2). Because the average time deviation among the three methods is small (~8.4%), we use the results from MCMCTREE-correlated-rate (its time estimates center within the three results) as our main dating results and illustrate them in Figure 2. Our molecular dating results suggest that the crown group of extant hynobiids originated in the middle Cretaceous, approximately 135 Ma (node a, Fig. 2). The major diversification of living hynobiids occurred in the late Eocene about 40 Ma (node c, Fig. 2). The

initial diversification of the most species-rich genus Hynobius is estimated ~ 17 Ma (node g, Fig. 2). The corresponding time estimates based on mitochondrial sequences for the two nodes is 62.5 Ma (Zhang et al., 2006) and 43 Ma (Li et al. 2011), about 50% and 135% older, respectively. Most divergence times estimated in the present study are close to those estimated with three nuclear loci (Zheng et al., 2011), but considerably younger than those from complete mitochondrial genomes (Zhang et al., 2006) (Fig. 2B). We agree with Zheng et al. (2011) on the reason for the discrepancy: because of substitutional saturation, fast evolving mitochondrial DNA sequences normally underestimate branch lengths for deep branches, making the time estimates systematically biased toward the ingroup root, which eventually produces considerably older dating results. However, we cannot simply think that dating of divergences using mitochondrial data is always less precise than using nuclear markers or that the dates derived from mitochondrial data are always overestimated. Mitochondrial substitutional saturation becomes an issue only for divergences occurring more than 15 million years ago (personal communication with Dr. Allan Larson). Therefore, when the range of divergence times being studied and the reference ingroup root are within this smaller timescale, for which mitochondrial substitutional saturation is minimal, divergence times based on mitochondrial data can still be precise. When the divergence times exceeding 15 Mya (herein hynobiid divergences can be as high as 140 million years ago) or the reference ingroup root is much older than the divergence times being estimated (very common in previous molecular dating

studies), for which mitochondrial substitutional saturation is severe, time estimates based on slowly evolving nuclear exons are more reliable because the nuclear exon data experienced much less substitutional saturation. However, using only a small number of nuclear genes (one to three, as in many recent studies) for molecular dating is not enough. Small sampling of nuclear genes usually causes sampling errors in collecting time signals and is also likely to produce imprecise time estimates. Actually, for most nodes in the timetree of hynobiids, our time estimates based on 29 nuclear genes are 10-20% older than previous results based on three nuclear genes (Zheng et al., 2011) and the 95% confidence intervals are also smaller (Table 2). To get more accurate timescales for many parts of vertebrate tree of life (especially for events occurring more than 15 million years ago), we recommend that dates estimated with a handful of mitochondrial genes or nuclear genes should be validated by further analyses using more nuclear loci.

Acknowledgments We thank Yingzhou Tian and Jianli Xiong for help in collecting samples XM1081 and XM1987. Xing Xing Shen provided valuable help in both lab works and data analyses. Many thanks to Editor Allan Larson for his thoughtful comments and help in improving the English usage of the early version of our manuscript. This work was supported by National Natural Science Foundation of China grants (31372172) and the National Science Fund for Excellent Young Scholars (No. 31322049) to P.Z.

Supplementary materials The sequence matrix and resulting tree files, and supplementary tables and figures associated with this study can be found in the online version of this article.

References AmphibiaWeb. 2014. Information on Amphibian Biology and Conservation. Available at http://amphibiaweb.org/. Accessed August 14, 2014. Ané C, Larget B, Baum D A, Smith SD, Rokas A. 2007. Bayesian estimation of concordance among gene trees. Mol. Biol. Evol. 24:412–426. Castresana J. 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17:540–552. dos Reis M, Zhu T, Yang Z. 2014. The impact of the rate prior on Bayesian estimation of divergence times with multiple Loci. Syst. Biol. 63:555–565. Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7:214. Edwards, J. L. 1976. Spinal nerves and their bearing on salamander phylogeny. J. Morphol. 148: 305–328. Edwards S V. 2009. Is a new and general theory of molecular systematics emerging? Evolution. 63:1–19. Gao K-Q, Shubin NH. 2003. Earliest known crown-group salamanders. 422:424–428. Guindon S et al. 2010. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst. Biol. 59:307–321.

Lanfear R, Calcott B, Ho SYW, Guindon S. 2012. Partitionfinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29:1695–1701. Lemmon EM, Lemmon AR. 2013. High-throughput genomic data in systematics and phylogenetics. Annu. Rev. Ecol. Evol. Syst. 44:99–121. Li J, Fu C, Lei G. 2011. Biogeographical consequences of Cenozoic tectonic events within East Asian margins: a case study of Hynobius biogeography. PLoS One. 6:e21506. Liu L, Yu L, Edwards S V. 2010. A maximum pseudo-likelihood approach for estimating species trees under the coalescent model. BMC Evol. Biol. 10:302. Löytynoja A, Goldman N. 2008. A model of evolution and structure for multiple sequence alignment. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 363:3913–3919. Near T.J., Eytan R.I., Dornburg A., Kuhn K.L., Moore J.A., Davis M.P., Wainwright P.C., Friedman M., Smith W.L. 2012. Resolution of ray-finned fish phylogeny and timing of diversification. Proc. Natl. Acad. Sci. U. S. A. 109: 13698-13703. Peng R et al. 2010. Rediscovery of Protohynobius puxiongensis (Caudata: Hynobiidae) and its phylogenetic position based on complete mitochondrial genomes. Mol. Phylogenet. Evol. 56:252–258. Pyron RA, Wiens JJ. 2011. A large-scale phylogeny of Amphibia including over 2800 species, and a revised classification of extant frogs, salamanders, and caecilians. Mol. Phylogenet. Evol. 61:543–583.

Regal, P. J. 1966. Feeding specializations and the classification of terrestrial salamanders. Evolution 20: 392–407. Roelants K., Gower D.J, Wilkinson M., Loader S.P., Biju S.D., Guillaume K., Moriau L., Bossuyt F. 2007. Global patterns of diversification in the history of modern amphibians. Proc. Natl. Acad. Sci. U. S. A. 104: 887-892. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Hohna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP 2012. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol 61: 539-542. Shen XX, Liang D, Feng YJ, Chen MY, Zhang P. 2013. A versatile and highly efficient toolkit including 102 nuclear markers for vertebrate phylogenomics, tested by resolving the higher level relationships of the caudata. Mol. Biol. Evol. 30:2235–2248. Shimodaira H, Hasegawa M. 2001. CONSEL: for assessing the confidence of phylogenetic tree selection. Bioinformatics. 17:1246–1247. Stamatakis A. 2006. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 22:2688–2690. Thorne JL, Kishino H. 2002. Divergence time and evolutionary rate estimation with multilocus data. Syst. Biol. 51:689–702.

Weisrock DW, Macey JR, Matsui M, Mulcahy DG 2013. Molecular phylogenetic reconstruction of the endemic Asian salamander family Hynobiidae (Amphibia, Caudata). Zootaxa 3626: 077–093. Zhang P, Chen YQ, Zhou H, Liu YF, Wang XL, Papenfuss TJ, Wake DB, Qu LH. 2006. Phylogeny, evolution, and biogeography of Asiatic salamanders (Hynobiidae). Proc. Natl. Acad. Sci. U. S. A. 103:7360–7365. Zheng Y, Peng R, Kuro-o M, Zeng X. 2011. Exploring patterns and extent of bias in estimating divergence time from mitochondrial DNA sequence data in a particular lineage: a case study of salamanders (order Caudata). Mol. Biol. Evol. 28:2521–2535.

Table Legends: Table 1. Summary information for the 29 nuclear genes used in this study.

Table 2. Detailed results of Bayesian molecular dating using MultiDivTime and MCMCTREE, and a comparison among shared nodes in previous studies. Letters for nodes correspond to Figure 2. Unit: one million years.

Figure Legends:

Figure 1. Phylogenetic relationships of hynobiid salamanders inferred from 29 nuclear loci (29,232 bp). The data set was analyzed under both concatenation framework (maximum-likelihood and Bayesian inference) and coalescent-based framework without data concatenation (MP-EST and BEAST). Nodes with bootstrap values > 90% and Bayesian posterior probabilities = 1.0 in the four different analyses are indicated as red filled squares. Otherwise, branches with letters have branch support values given beside the tree. Tree topology and branch lengths were from the concatenated ML analysis. Outgroup species are not shown.

Figure 2. (A) Time-calibrated phylogeny of hynobiids inferred from 29 nuclear genes. The chronogram was based on the results estimated by the method of MCMCTREE (correlated rates). The light-blue bars through the nodes indicate 95% credibility intervals. Detailed time estimates for nodes with letter labels are given in Table 2. (B) Comparison of divergence-time estimates for eight nodes shared across three studies. The circle within boxes represents the mean of the posterior estimate and the whiskers mark the upper and lower 95% highest posterior density of the age estimates. The comparison shows that our new time estimates are largely congruent with previous results based on three nuclear genes, but considerably younger than those estimated by complete mitochondrial genomes.

Onychodactylus fischeri O. zhangyapingi Paradactylodon mustersi Ranodon sibiricus Pachyhynobius shangchengensis Hynobius retardatus H. kimurae H. amjiensis H. chinensis

A

H. guabangshanensis H. leechii H. nigrescens

B C Node

H. tokyoensis Salamandrella keyserlingii

Support value Concatenated

Concatenated

Species-tree

Species-tree

ML

BI

MP-EST

BEAST

A

100%

1.0

93%

0.99

B

81%

1.0

57%

0.87

C

66%

0.98

85%

0.91

D

94%

1.0

64%

1.0

E

59%

0.96

-

-

F

99%

1.0

54%

1.0

G

59%

0.96

77%

0.92

Batrachuperus yenyuanensis B. pinchonii D

B. longdongensis E F

B. karlschmidti B. tibetanus Liua tsinpaensis L. shihi

G

Pseudohynobius puxiongensis P. flavomaculatus

0.01

P. shuichengensis

Jurassic

Cretaceous

160

140

(B)

160

120

Cenozoic

100

80

60

40

20

0 Liua shihi

w a

140

Liua tsinpaensis

t

mtDNA genome (Zhang et al., 2006) 3 nuclear loci (Zheng et al., 2011) 29 nuclear loci (this study)

Pseudohynobius flavomaculatus

v u

Pseudohynobius shuichengensis

Divergence time (Ma)

120

Pseudohynobius puxiongensis

o 100

80

c 60

e d

s

Batrachuperus karlschmidti

r

Batrachuperus tibetanus

q

n

Batrachuperus pinchonii

p

f

Batrachuperus longdongensis

o

40

Batrachuperus yenyuanensis p

Salamandrella keyserlingii

w

20

f 0

m

Hynobius tokyoensis

l

Hynobius nigrescens Hynobius leechii

i

Hynobius chinensis

k h

e

Hynobius guabangshanensis j

g

Hynobius amjiensis Hynobius kimurae

c

(A)

Hynobius retardatus Pachyhynobius shangchengensis

a

d

Ranodon sibiricus Paradactylodon mustersi

Root

Onychodactylus zhangyapingi

b

Onychodactylus fischeri Andrias davidianus

160

140

120

100

80

60

40

20

0

Table 1. Summary information for the 29 nuclear genes used in this study. Locus

Alignment length(bp)

Number of taxa sequenced

Number of variable sites

Overall mean pairwise P distance

Best-fitted model

PTF CAND1 DET1

549 1155 699

31 31 31

140 361 223

0.046 0.060 0.066

HKY+G HKY+G K80+G

DISP1 DISP2 DNAH3 DSEL ENC1 EXTL3 FAT1 FAT2

954 942 915 1263 1083 1215 1503 936

31 31 31 31 31 31 30 30

252 356 346 479 319 386 653 420

0.045 0.075 0.07 0.073 0.066 0.058 0.080 0.083

HKY+G K80+G HKY+G HKY+G K80+I+G K80+I+G HKY+G GTR+G

FAT4 HYP KBTBD2 KCNF1 KIAA1239 LIG4 LPHN2 MB21D2

711 1260 1071 762 1275 1005 579 1008

31 25 31 31 31 31 31 31

276 441 285 269 383 365 159 316

0.080 0.077 0.050 0.082 0.058 0.067 0.052 0.071

K80+G SYM+G HKY+G K80+G HKY+G HKY+G K80+G K80+G

PANX2 PDP1 PIK3CG PPL RAG1 RAG2 SACS STLITRK1

696 1038 927 1239 1275 873 1080 1143

31 31 31 31 31 31 30 30

226 308 381 500 465 417 310 407

0.072 0.060 0.085 0.080 0.077 0.099 0.059 0.074

K80+G HKY+I+G K80+G GTR+G K80+I+G K80+G HKY+G GTR+G

ZBED4 ZHX2

1002 1074

31 31

264 474

0.045 0.085

HKY+G HKY+I

Table 2. Detail results of Bayesian molecular dating using MultiDivTime and MCMCTREE and a comparison among shared node in previous studies. Letters for nodes are corresponding to Figure 2. Unit: one million year.

MCMCTREE

MCMCTREE

(Independent rates)

(correlated rates)

mean (95% CI)

mean (95% CI)

mean (95% CI)

Root*

152.0 (145.2-163.5)

153.7 (144.9-164.7)

157.1 (145.6-165.3)

a

121.5 (112.3-132.9)

126.8 (113.4-141.8)

135.1 (120.2-150.3) 117 (92-141)

MultiDivTime Node

Zheng et al. 2011

Zhang et al. 2006

mean (95% CI)

mean (95% CI)

110.7 (106.4-114.9)

b

16.2 (11.5-21.7)

5.8 (4.1-7.7)

7.5 (5.1-10.1)

c

34.1 (29.5-39.4)

41.0 (36.1-46.5)

40.2 (34.5-46.2)

36 (27-45)

62.5 (59.7-65.5)

d

22.9 (19.3-27.0)

28.0 (23.8-32.7)

27.0 (22.7-31.9)

21 (13-29)

48.3 (45.7-50.9)

e

30.0 (25.8-34.8)

36.8 (32.4-41.8)

36.2 (31.1-41.7)

30 (22-38)

57.8 (55.2-60.5)

f

28.6 (24.5-33.2)

35.6 (31.3-40.4)

35.1 (30.1-40.5)

22 (16-28)

52.5 (50.0-54.9)

g

14.1 (11.6-16.9)

17.5 (14.9-20.4)

16.6 (13.9-19.7)

14 (10-18)

h

13.4 (11.0-16.1)

16.9 (14.4-19.7)

16.2 (13.5-19.2)

i

9.6 (7.6-11.8)

11.1 (9.3-13.1)

10.5 (8.6-12.6)

24.8 (23.3-26.5)

j

4.1 (2.8-5.6)

3.6 (2.5-4.8)

3.5 (2.4-4.7)

19.4 (16.6-20.8)

k

3.4 (2.2-4.8)

3.4 (2.3-4.6)

3.3 (2.2-4.5)

l

8.8 (7.0-10.9)

10.8 (9.0-12.8)

10.2 (8.3-12.3)

m

8.0 (6.3-10.0)

10.4 (8.7-12.4)

9.9 (8.0-12.0)

n

26.7 (22.7-31.1)

33.8 (29.6-38.5)

33.3 (28.4-38.6)

o

15.9 (13.1-19.0)

15.4 (13.0-18.1)

16.2 (13.5-19.3)

14 (10-19)

42.9 (40.2-44.6)

p

8.4 (6.4-10.6)

7.1 (5.4-8.9)

7.3 (5.5-9.4)

10 (6-14)

24.3 (22.6-26.0)

q

6.1 (4.5-8.0)

4.1 (2.9-5.5)

4.5 (3.1-6.0)

r

4.6 (3.1-6.2)

3.8 (2.6-5.1)

4.0 (2.7-5.5)

s

2.9 (1.5-4.3)

3.1 (1.9-4.4)

3.1 (1.7-4.5)

48.5 (46.1-51.0)

t

14.0 (11.4-17.0)

15.0 (12.7-17.6)

15.8 (13.1-18.8)

12 (8-16)

u

6.9 (5.1-8.8)

6.2 (4.5-8.0)

6.7 (5.0-8.7)

5 (3-8)

v

4.7 (3.2-6.4)

4.6 (3.2-6.2)

5.0 (3.5-6.7)

3 (1-5)

w

5.3 (3.7-7.0)

4.2 (2.7-5.8)

4.6 (3.2-6.2)

5 (2-8)

* Calibration points.

21.5 (19.8-23.1)


29 nuclear genes (a total number of 29,232 bp) were used in this study
Concatenation and specie-tree analyses yielded congruent and robust results
The phylogenetic positions for Pachyhynobius and Salamandrella are determined
Previous time estimates for hynobiids based on mtDNAs are overestimated

Previous studies (conflicting results)

Salamandrella Pachyhynobius

Pachyhynobius

Ranodon Ranodon

Pachyhynobius

Hynobius

Pseudohynobius Salamandrella

Salamandrella

Paradactylodon

Paradactylodon

Liua

Salamandrella

Onychodactylus

Onychodactylus

Batrachuperus

Pachyhynobius

This study (29 nuclear genes)

Hynobius

100/1/99/1 100/1/97/1 94/1/64/1 100/1/100/1

Salamandrella Batrachuperus Liua Pseudohynobius