Comparative phylogeography of bamboo bats of the genus Tylonycteris ( Chiroptera , Vespertilionidae ) in Southeast Asia

1,5,6 Institute of Ecology and Biological Resources, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet Road, Cau Giay District, Hanoi, Vietnam. 1,7,8 Muséum national d’Histoire naturelle, Service de Systématique Moléculaire, UMS 2700, CP 26, 43, Rue Cuvier, 75005 Paris, France. 1,8 Institut de Systématique, Evolution, Biodiversité (ISYEB), Sorbonne Universités, UMR 7205 MNHN CNRS UPMC, Muséum national d'Histoire naturelle, CP 51, 55, Rue Buffon, 75005 Paris, France. 2 Department of Zoology, Hungarian Natural History Museum, Baross u. 13, 1088 Budapest, Hungary. 3 Department of Mammalogy and Ornithology, Natural History Museum of Geneva, Route de Malagnou 1, BP 6434, 1211 Geneva 6, Switzerland. 4 Fauna & Flora International, Cambodia Programme, 19 Street 360, BKK 1, Chamkarmorn, Phnom Penh, Cambodia.


Introduction
Bamboo bats of the genus Tylonycteris Peters, 1872 (Chiroptera, Vespertilionidae) are small-sized bats (weight: 3-10 g; forearm length: 22-32 mm) characterized by a dorsoventrally flattened skull and well-developed, fleshy pads at the base of the thumb and on the sole of the foot (Tate 1942).These morphological features are thought to be adaptations for roosting in small cavities with smooth surfaces such as the internodes of bamboo stalks or narrow crevices in trees and rocks (Feng et al. 2008;Medway & Marshall 1970;Thewissen & Etnier 1995).Classically, the genus was regarded as containing only two species, T. pachypus (Temminck, 1840) and T. robustula Thomas, 1915, with several further taxa included as subspecies (Fig. 1) (Tate 1942;Simmon 2005).More recently, Feng et al. (2008) described a third species, T. pygmaea Feng, Li & Wang, 2008, which is smaller than its congeners.Whereas T. pygmaea is considered to be endemic to Yunnan Province in southern China, the two other species have much more extensive geographic ranges that greatly overlap in Southeast Asia (Fig. 1).The range of T. pachypus is even more extended in the north and west, with apparently two isolated populations recorded in the Chinese provinces of Sichuan and Chongqing, and in southern India around the Western Ghats (Fig. 1) (Bates et al. 2008a(Bates et al. , 2008b)).
Bats of the genus Tylonycteris are usually associated with bamboo groves in both intact and disturbed habitats at elevations ranging from lowland up to 1500 m (Bates et al. 2008a(Bates et al. , 2008b;;Medway & Marshall 1970).All species can often be found roosting within bamboo internodes in colonies of up to 40 individuals, entering through narrow vertical slits created by long-horned beetles (Medway & Marshall 1970;Zhang et al. 2007).Females are gregarious, whereas males tend to be more solitary (Medway & Marshall 1972).At least occasionally, different species of Tylonycteris can occupy the same bamboo chamber (Feng et al. 2008;Medway & Marshall 1970).Zhang et al. (2005) determined that the diet of Tylonycteris bats is mainly composed of insects of the order Hymenoptera, which were regarded as significant pests of bamboo.Interestingly, the highest richness in bamboo species in the Asia-Pacific Region has been recorded in southern China (Bystriakova et al. 2003a;Yuming et al. 2004), the only region where the three species of Tylonycteris occur in sympatry.Taken together, these aspects indicate that Tylonycteris species are strongly dependent on bamboo vegetation.Thomas, 1915(Bates et al. 2008a, 2008b).C. Localities of Tylonycteris specimens included in this study.Triangles and circles refer to T. pachypus and T. robustula, respectively; the colours indicate the mtDNA haplogroups found in the Bayesian analyses of COI and Cytb sequences (see Fig. 2

for details).
European Journal of Taxonomy 274: 1-38 (2017) 4 T. p. pachypus found in Peninsular Malaysia and Indonesia, suggesting that T. pachypus might be a species complex.
We present here a comprehensive comparative phylogeographic investigation of Tylonycteris spp.based on molecular and morphological data.Two mitochondrial markers and seven nuclear markers were sequenced to address the following questions: (1) how many mtDNA haplogroups exist in mainland Southeast Asia; (2) are these haplogroups corroborated by nuclear DNA markers; and (3) do these genetically defined groups correspond to distinct morphotypes?The ultimate goal of this study was to improve our understanding of the biogeography and taxonomy of the genus Tylonycteris.

Data sampling for genetic analyses
For this study, 63 tissue samples were collected during field expeditions within Indochina and from historical specimens housed in museum collections.Vouchers were deposited in the following institutions: the Institute of Ecology and Biological Resources (IEBR, Hanoi, Vietnam), the Muséum national d'Histoire naturelle (MNHN, Paris, France), the Muséum d'histoire naturelle de Genève (MHNG, Geneva, Switzerland), the Hungarian Natural History Museum (HNHM, Budapest, Hungary), the Centre for Biodiversity Conservation (CBC, Royal University of Phnom Penh, Cambodia), Rijksmuseum van Natuurlijke Histoire (RMNH, Leiden, Netherlands; now Naturalis Biodiversity Center) and the Zoological Museum Amsterdam (ZMA, Amsterdam, Netherlands).In addition, we included DNA sequences of Tylonycteris available in the GenBank nucleotide database.The origins of all samples are represented in Fig. 1 and detailed in Appendices 1 and 2. Bats were captured using mist nets (Ecotone, Poland) and four-bank harp traps.They were measured, photographed and initially identified following morphological criteria provided by Bates & Harrison (1997), Borisenko & Kruskop (2003) and Francis (2008).Tissue samples for genetic study were collected from chest muscles of voucher specimens or from the patagium (biopsy punches; 3 mm diameter) of released bats.Tissue samples were preserved in 95% ethanol and stored at -20°C.

DNA extraction, amplification, sequencing
Total DNA was extracted from muscle or patagium samples using the QIAGEN DNAeasy Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol.Two mitochondrial genes were sequenced for this study: the 5' fragment of the COI gene and the complete cytochrome b (Cytb) gene.The primers used for PCR amplification of mitochondrial genes were UTyrLA and C1L705 for COI (Hassanin et al. 2012), and Thr-CH (Hassanin et al. 2014) and the newly designed Glu-CH (5'-AAY CAC CGT TGT AYT TCA ACT A-3') for Cytb.From DNA extracted from five old museum vouchers, it was not possible to amplify long PCR fragments (> 300 nt).Therefore, two further primer sets were used to amplify and sequence two overlapping fragments of COI: UTyrLA and Tyl-L122 (5'-AGY ART GCT CCT GGY TGA CC-3'), and Tyl-U86 (5'-GGT GCC TGA GCT GGC ATA GT-3') and Tyl-L266 (5'-TCG GGG RAA TGC TAT ATC AG-3').
The following seven nuclear markers were also sequenced for three outgroup species and 14 samples of Tylonycteris representing the divergent mtDNA haplotypes of Indochina: intron 2 of CHPF2 (chondroitin polymerizing factor 2), intron 6 of HDAC1 (histone deacetylase 1), intron 10 of HDAC2 (histone deacetylase 2), intron 2 of PABPN1 (poly(A) binding protein, nuclear 1), intron 6 of RIOK3 (RIO kinase 3), intron 9 of TUFM (elongation factor Tu, mitochondrial precursor), and intron 6 of ZFYVE27 (zinc finger, FYVE domain containing 27).The primers used for PCR amplifications of nuclear introns are detailed in Hassanin et al. (2013).For Sumatran specimens, which were collected more than one hundred years ago, it was not possible to obtain successful PCR amplification of nuclear markers.Amplifications were done in 20 µl using 3 µl of Buffer 10X with MgCl 2 , 2 µl of dNTP (6.6 mM), 0.12 µl of Taq DNA polymerase (2.5 U, Qiagen, Hilden, Germany) and 0.5-1.0µl of the two primers at 10 µM.The standard PCR conditions were as follows: 4 min at 95°C; 5 cycles of denaturation/ annealing/extension for 45 s at 95°C, 1 min at 60°C and 1 min at 72°C, followed by 30 cycles of 30 s at 95°C, 45 s at 55°C, and 1 min at 72°C, followed by 10 min at 72°C.PCR products were resolved by electrophoresis on a 1.5% agarose gel stained with ethidium bromide and visualized under UV light.Both strands of PCR products were sequenced using Sanger sequencing on an ABI 3730 automatic sequencer at the Centre National de Séquençage (Genoscope) in Evry (France).The sequences were edited and assembled using Codoncode Alignment v. 3.7.1 (CodonCode Corporation) and Sequencher v. 5.0 (Gene Codes Corporation).Heterozygous positions (double peaks) were scored using the IUPAC ambiguity codes.Sequences generated for this study were deposited in the GenBank database (accession numbers KX496340-KX496537; Appendix 1).

Phylogeographic analyses based on mitochondrial sequences
The 63 COI and 19 Cytb sequences newly generated in this study were compared to the 38 COI and four Cytb sequences available for Tylonycteris in GenBank.DNA sequences were aligned with Se-Al v. 2.0a11 (Rambaut 2002).Phylogenetic analyses of bats of the genus Tylonycteris were initially performed using two separate datasets: ( 1) COI (101 taxa and 728 nt), and (2) Cytb (23 taxa and 1140 nt).Three newly sequenced species, representing three different genera of the subfamily Vespertilioninae (Pipistrellus cf.javanicus Gray, 1838, Eptesicus sp. and Hypsugo pulveratus (Peters, 1870)), were chosen as outgroups on the basis of previous studies (Jones et al. 2002;Roehrs et al. 2010).
The Bayesian approach was used to reconstruct phylogenetic relationships.The best-fitting model of sequence evolution was selected under jModelTest (Posada 2008) using the Akaike Information Criterion (AIC).Bayesian analyses were then conducted with MrBayes v. 3.2.1 (Ronquist et al. 2012) using the selected HKY + G and GTR + I + G models, respectively for COI and Cytb datasets.The posterior probabilities (PP) were calculated using four independent Markov chains run for 10 000 000 Metropolis-coupled MCMC generations, with trees sampled every 1000 generations and a burn-in of 25%.Pairwise genetic distances between groups were calculated with PAUP* v. 4b10 (Swofford 2003) using the Kimura 2-parameter (K2P) correction.

Phylogenetic analyses based on eight independent markers
The phylogeny of Tylonycteris was further investigated using a reduced sample of 14 specimens (plus three outgroups) sequenced for multiple loci, including seven nuclear introns and a combination of two mtDNA markers (Appendix 1).These analyses were performed to test possible discordance between the phylogenetic signals extracted from independent markers.We did not obtain high quality sequences for three PCR products: CHPF2 for two specimens of T. pachypus (T5009 and VN11-1138) and PABPNI for Eptesicus sp.(VN11-0076) (Appendix 1).These poor-quality sequences were not included in the alignments and coded as missing data (Ns) in the multilocus analyses.Accordingly, 10 datasets were analyzed: supermatrix (combining all the nine genes; 7572 nt), nuDNA (combining all the seven nuclear introns; 5604 nt), mtDNA (COI + Cytb; 1868 nt), CHPF2 (858 nt), HDAC1 (1128 nt), HDAC2 (639 nt), PABPNI (677 nt), RIOK3 (915 nt), TUFM (655 nt) and ZFYVE27 (732 nt).DNA alignments were done with Se-Al v. 2.0a11 (Rambaut 2002).A few gaps were included in the alignments of the nuclear introns, but their positions were not found to be ambiguous.The models of nucleotide evolution were selected under jModeltest (Posada 2008): the GTR + I + G model for mtDNA markers, the HYK + G model for CHPF2, PABPNI and ZFYVE27, and the HKY model for HDAC1, HDAC2, RIOK3 and TUFM.
A partitioned Bayesian approach was used to account for the combination of markers with contrasted molecular properties.The mtDNA dataset was run using a GTR + I + G model for each of the three codon positions of the two mitochondrial genes; the concatenation of seven nuclear introns and the supermatrix were run using the selected model for each gene partition.The indels shared by at least two taxa and without ambiguity in the position of the gaps were coded as additional characters ("1": insertion; "0": deletion) and analyzed as a separate partition in the Bayesian analyses.The posterior probabilities (PP) were calculated using four independent Markov chains run for 10 000 000 Metropolis-coupled MCMC generations, with tree sampling every 1000 generations, and a burn-in of 25%.
Bootstrap percentages (BP) were computed by PAUP* v. 4b10 (Swofford 2003) after 1000 replicates, using the GTR + I + G model selected by jModelTest for the supermatrix dataset.The results obtained from the separate Bayesian analyses of the eight independent molecular markers (mtDNA and the seven nuclear introns) were also analyzed for congruence using the SuperTRI method (Ropiquet et al. 2009).The lists of bipartitions obtained from the eight Bayesian analyses were transformed into a weighted binary matrix for supertree construction using SuperTRI v. 57 (available from http://www.normalesup.org/~bli/Programs/programs.html).Each binary character corresponds to a node, which was weighted according to its frequency of occurrence in one of the eight lists of bipartitions.In that way, the SuperTRI method takes into account both principal and secondary signals, because all phylogenetic hypotheses found during the eight separate analyses are represented in the weighted binary matrix used for supertree construction.The reliability of the nodes was assessed using three different measures.The first value is the Supertree Bootstrap Percentage (SBP), which was calculated under PAUP* v. 4b10 after 1000 BP replicates of the weighted binary matrix reconstructed with SuperTRI (941 characters; heuristic search).The second value is the ''Mean Posterior Probability'' (MPP) calculated using the lists of bipartitions obtained from Bayesian analyses of the eight datasets.The third value is the index of reproducibility (Rep), which is the ratio of the number of datasets supporting the node of interest to the total number of datasets.The MPP and Rep values were directly calculated on SuperTRI v. 57.All SuperTRI values were reported on the Bayesian tree obtained from the supermatrix analysis.

Molecular Dating
Divergence times were estimated with the Bayesian approach implemented in BEAST v. 2. 1.3 (Drummond et al. 2012) using either a COI alignment of 33 taxa and 291 nt or a Cytb alignment of 19 taxa and 1140 nt.As no calibration point (fossil record or biogeographic event) is available for Tylonycteris, we employed a range of nucleotide substitution rates used for mammals (Arbogast & Slowinski 1998;Mao et al. 2010).We used a first mutation rate (R1) of 0.015 per site per lineage per Myr with a lower boundary of 0.01 and an upper boundary of 0.02, a second mutation rate (R2) of 0.02 ± 0.005 per site per lineage per Myr and a third mutation rate (R3) of 0.025 ± 0.005 per site per lineage per Myr.We tested R1 and R2 for the CO1 dataset, and R2 and R3 for Cytb dataset.We applied a GTR + I + G model of evolution (based on jModelTest) and a relaxed-clock model with an uncorrelated log normal distribution for substitution rate.Node ages were estimated using a Yule speciation prior and 10 8 generations, with tree sampling every 1000 generations and a burn-in of 25 %.Adequacy of chain mixing and MCMC chain convergence were assessed using the ESS values in Tracer v. 1.6 (available in the BEAST package).The chronograms were generated with TreeAnnotator v. 1.8.2 (also available in the BEAST package) and visualized with FigTree v. 1.4.1 (Rambaut 2009).

Morphological analysis
For the morphological analysis, 62 adult specimens of Tylonycteris spp.were examined (Appendix 1).Besides mass (expressed in grams), the following external measurements were taken from living bats or museum specimens to the nearest 0.1 mm: forearm length (FA) -from the elbow to the wrist with both joints folded; head and body length (HB) -from the tip of the face/chin to the anus; tail length (Tail) -from the anus to the tip of the tail; ear length (Ear) -from the base of the ear, where it attaches to the head, to the tip of the pinna.Craniodental measurements were taken to the nearest 0.01 mm with the use of digital calipers under a stereo microscope and included: GLS -greatest length of skull, from the anterior rim of the alveolus of the 1 st upper incisor to the most posteriorly projecting point of the occipital region; CCL -condylo-canine length, from the exoccipital condyle to the most anterior part of the canine; UCI -from the anterior rim of the alveolus of the first upper incisor to the posterior rim of the alveolus of the canine; CC -greatest width across the upper canines from their labial borders; M 3 M 3 -greatest width across the crowns of the last upper molars from their labial borders; IC -least width of the interorbital constriction; MB -greatest distance across the mastoid region; BW -greatest width of the braincase; CM 3 -maxillary toothrow length, from the anterior of the upper canine to the posterior of the crown of the 3 rd molar; ML -mandible length, from the anterior rim of the alveolus of the 1 st lower incisor to the most posterior part of the condyle; and CM 3 -mandibular toothrow length, from the anterior of the lower canine to the posterior of the crown of the 3 rd lower molar.

Data analyses
Bats were classified into six geographic groups on the basis of our molecular results.Two distinct principle component analyses (PCAs) were done on the craniodental characters with PAST (Hammer et al. 2001) and irrespective of sexes: (1) the log-transformed raw measurements to assess an overall size factor (usually PC1; Barlow et al. 1997;Lindenfors et al. 2007) and ( 2) log-transformed standardized data (raw score/geometric mean) to assess shape factors (Jungers et al. 1995).The statistically significant difference in mean values of craniodental measurements and PC mean scores between different groups were then tested by one-way analysis of variance (ANOVA, p ≤ 0.05) followed by Turkey HSD post-hoc test for unequal sample sizes (Zar 1999).

Results derived from the analysis of COI sequences
We sequenced a COI fragment of 728 nt for 55 individuals of Tylonycteris and three outgroup taxa collected across Indochina and adjacent areas.Five additional COI sequences of 291 nt were also obtained from old specimens collected from Sumatra (Appendix 1).Our sequences were aligned with all other COI sequences available in the GenBank database, producing a final alignment containing 101 specimens.The Bayesian analyses reconstructed from the two COI alignments of either 728 nt (Fig. 2) or 291 nt (Appendix 3) showed very similar topologies.Both analyses support the monophyly of the genus Tylonycteris (PP = 1).As noted by Huang et al. (2014), one specimen labelled as T. robustula in GenBank (accession number HM914921) appears within the clade containing all sequences of T. pachypus.Examination of the voucher, ZMMU S-186637 held at the Zoological Museum of Moscow State University, indicates that its cranium clearly fits the diagnosis of T. pachypus (S.Kruskop, in litt.) and that it should be relabelled accordingly.Otherwise, our results strongly support the monophyly of both species of Tylonycteris.In the COI trees (Fig. 2; Appendix 3), the T. pachypus clade is further subdivided into three divergent and highly supported geographic subclades (PP = 1) named Tp1 (Sumatra), Tp2 (southern Indochina) and Tp3 (northern Indochina).A similar phylogeographic pattern is also detected in T. robustula, with three strongly supported geographic subclades named Tr1 (Sumatra), Tr2 (southern Indochina + northwestern India + Peninsular Malaysia) and Tr3 (northern Indochina).Within Indochina, the area of separation between the northern and southern Indochinese subclades corresponds to the region located between Vu Quang (Ha Tinh Province, Vietnam) and Ban Houana (Khammouane Province, Laos) for T. pachypus, whereas it corresponds to the region situated between Nam Et (Houaphan Province, Laos) and Vu Quang (Vietnam) for T. robustula (Fig. 1).

Results derived from the analysis of Cytb sequences
The Bayesian analysis of the newly generated Cytb sequences (n = 19) and those of Tylonycteris downloaded from GenBank (n = 4) produced a topology (Fig. 2) similar to those obtained with the COI gene tree, with two exceptions: (1) T. robustula is paraphyletic because one specimen collected Fig. 2. Bayesian tree obtained from the analyses of COI and Cytb genes.The values on the nodes indicate posterior probabilities (PP).In the COI tree, PPs were calculated from two DNA alignments of CO1 sequences (728 nt or 291 nt; see text for details).The sequence HM914921 was obtained from a specimen originally identified as T. robustula.In the Cytb tree, the sequence EU521635 was obtained from a specimen originally identified as T. robustula.(2) two sequences of T. pachypus from Guangdong and Guangxi (China) are nested within the northern Indochinese subclade (Tp3) of T. pachypus (PP = 0.7) as expected, but a sequence of T. pachypus from Hong Kong appears as sister to the southern Indochina subclade (Tp2) (PP = 1) (Fig. 2).These results are consistent with the phylogeographic pattern reported in Huang et al. (2014).

Genetic distances
The comparisons of interspecific genetic variations, as estimated by K2P distances, indicate that the specimens from Sumatra differ from those of mainland Southeast Asia by 6.0-6.1% and 5.7-7.5% in partial COI sequences (calculated from an alignment of only 291 nt) for T. pachypus and T. robustula, respectively (Appendix 5).The K2P distances for Cytb sequences (404 nt) between the specimen of T. robustula from Borneo and those representing the Tp2 and Tp3 subclades from Southeast Asia's mainland (Indochina and southern China including Guangdong, Guangxi, and Hong Kong) range from 5.7 to 6.4%.This is significantly smaller than the distance of 12.4-13.8%measured for other specimens assigned to T. robustula (Appendix 5).Within Indochina, the mean K2P distances inferred from the alignments of COI (728 nt) and Cytb (1140 nt) sequences are 2.8% (2.2-3.5%) and 2.8% (2.6-3.0%) between Tp2 and Tp3 of T. pachypus, and 6.5% (5.5-7.4%) and 9.5% (8.6-10.4%) between Tr2 and Tr3 of T. robustula, respectively.The maximal intraspecific distances in Cytb sequences between the specimen from Hong Kong and those from Tr2 and Tr3 are 2.0% and 3.2% (3.1-3.2%),respectively (data not shown in Appendix 5).

Supermatrix and SuperTRI analyses of eight independent DNA markers
The Bayesian trees reconstructed from the concatenation of the seven nuDNA introns (5604 nt and 54 indels) or based on the supermatrix of 7526 characters (mtDNA + nuDNA; 7472 nt and 54 indels) are shown in Fig. 3 and Appendix 4, respectively.In the tree of Fig. 3, we report the bootstrap proportions obtained from the ML analyses (Appendix 6), as well as the results of the SuperTRI analyses (SBP, MPP, and Rep, Appendix 7).All these analyses support with maximal values of robustness the monophyly of the genus Tylonycteris and that of the two species, T. pachypus and T. robustula (Fig. 3).Both species can be diagnosed by several indels in the nuclear dataset: 10 for T. pachypus and three for T. robustula (Fig. 3).For instance, all specimens of T. pachypus share a large insertion of ca 250 nt in RIOK3 and a deletion of 47 nt in HDAC2, while a deletion of two bases (CA) in CHPF2 was found in all specimens of T. robustula (Fig. 3).
Within T. robustula, the two geographic mtDNA haplogroups corresponding to Tr2 and Tr3 in Indochina were recovered with high support values (PP nuDNA / PP supermatrix = 1; BP = 100; SBP > 95; MPP > 0.54; Rep > 0.5; Fig. 3).Other relationships within T. robustula are not robust.Within T. pachypus, no substructure was consistent.The Tp2 subclade was recovered in the supermatrix analyses with modest support and SuperTRI analyses revealed that this node was only supported by the mtDNA dataset.The Tp3 subclade was not found to be monophyletic in supermatrix and superTRI analyses, but the weak support values indicate a lack of resolution rather than real discordance between the datasets.
Interspecific divergences, as estimated from K2P distances, were calculated using the concatenation of the seven nuclear genes (5604 nt).The distances between T. pachypus and T. robustula range from 1.6 to 1.8%.The divergences between specimens of northern and southern Indochina range from 0.4 to 0.6% for T. robustula and from 0 to 0.2% for T. pachypus.The maximum intraspecific variation within each geographical population is 0.2% for both species (Appendix 5).

Molecular divergence time estimates
Our molecular dating estimations based on two different mtDNA alignments (COI and Cytb) and three different rates of substitution (Table 1) suggest that the common ancestor of extant Tylonycteris spp. was already present during the Pliocene epoch (6.56-3.84Mya), and that the diversification of major geographic clades within both T. pachypus and T. robustula occurred during the Early Pleistocene.

Morphological comparisons
Most craniodental measurements of available T. pachypus specimens are significantly smaller than those of T. robustula (ANOVA, p ≤ 0.05); whereas those of the three molecular groups within each taxon are overlapping (Table 2).This is reflected in univariate variation of GLS measurements (Fig. 4A) or on the multivariate component of PC*1 based on raw data.This component accounted for 86.9% of total variance and was correlated positively with all characters (Table 3), confirming that it represents an overall size factor (Fig. 4B).The range of PC*1 scores of T. pachypus (between -2.0 and -0.58) is significantly smaller than that of T. robustula (between -0.38 and 1.33) (ANOVA, p ≤ 0.05).
The PCA of log-transformed standardized data, which better reflects shape factors, reveals that the two first PCs, PC1 and PC2, accounted for 37.2 and 19.3% of total variance, respectively (Table 3), with a Fig. 3. Bayesian tree reconstructed from the analysis of the seven concatenated nuDNA introns.The seven independent markers are CHPF2 (858 nt), HDAC1 (1128 nt), HDAC2 (639 nt), PABPNI (677 nt), RIOK3 (915 nt), TUFM (655 nt) and ZFYVE27 (732 nt).The values indicated on the branches are the Bayesian posterior probabilities (PP nuDNA and PP supermatrix ), maximum likelihood bootstrap (BP), Supertree Bootstrap percentage (SBP), Mean posterior probability (MPP) and Reproducibility index (Rep).An asterisk indicates that the node was supported by maximal values of robustness (PP nuDNA = PP supermatrix = 1; BP = 100).The letter ''X'' indicates that the node was not found in the analysis.The position and nature of all diagnostic indels (i: insertion; d: deletion) shared by at least two taxa in the alignments of nuclear genes are indicated in boxes.

Tr2
Tr3 significant difference in shape components between Tylonycteris taxa (ANOVA, p ≤ 0.05).In the plot of PC1 against PC2 (Fig. 4C), we found an overlap between T. pachypus and T. robustula.However, this multivariate projection also shows that, within T. pachypus, bats of Tp1 form a cluster clearly separated from bats of Tp2 and Tp3.Within T. robustula, morphological overlap between haplogroups is more extensive, but most Tr2 individuals tend to separate from Tr1 and Tr3 ones.

Systematics of the genus Tylonycteris
In earlier taxonomic accounts (e.g., Tate 1942 or Medway 1973), the three distinct species T. fulvida (sometimes incorrectly referred to as rubidus) (type locality [t.l.]: Schwe Gyen, Myanmar), T. meyeri Peters, 1872 (t.l.: South Luzon) and T. aurex Thomas, 1951 (t.l.: Astoli, south of Mumbai, India) were treated as subspecies of T. pachypus (t.l.: Java), while T. malayana Chasen, 1940 (t.l.: Peninsular Malaysia) was considered as a subspecies of T. robustula (t.l.: Borneo).Most subsequent authors followed these taxonomic recommendations (Corbet & Hill 1992;Simmons 2005;Bates et al. 2008aBates et al. , 2008b;;Francis 2008).Previous karyological and genetic data (Francis et al. 2010;Huang et al. 2011) and our new phylogenetic analyses, however, support a clear division between populations of Tylonycteris spp.from mainland Southeast Asia vs those of the Sundaic islands, such as Sumatra and Borneo.We propose therefore to restrict the species names T. pachypus and T. robustula to populations from the Sundaland, where the type specimens were collected (Fig. 1), and to revalidate the following two species names previously used for populations occurring in mainland Southeast Asia (Tate 1942;Medway 1973): the first is T. fulvida, which was originally described as Scotophilus fulvidus from Sittang River (Myanmar) and should be applied to continental bats of the T. pachypus complex; the second is T. malayana, which was originally described from Perak (Peninsular Malaysia) and should be applied to populations of the T. robustula complex found in Peninsular Malaysia, southern Indochina, and northern India.In contrast, all bats of the T. robustula species complex endemic to northern Indochina should be included in the new species described below.Representatives from other isolated populations, notably from the Philippines and the Western Ghats in India would need further scrutiny to assess their taxonomic rank, but certainly represent further distinct taxa.Tylonycteris robustula Thomas, 1915 (partim): 227.

Etymology
The specific epithet refers to the current restricted occurrence of the new species in north-eastern Laos and northern Vietnam (Fig. 1).The Vietnamese portion of this region was previously called "Tonkin"

Referred material
Specimens identified as T. robustula collected from Na Don, Phuong Vien, Cho Don (Bac Kan Province) and Na Hang Nature Reserve (Tuyen Quang Province) (Appendix 1) are also referred to T. tonkinensis sp.nov.

Description
A member of the T. robustula species complex comprising representatives of the Tr3 haplogroup found in northern Indochina.Externally, individuals are small, with a forearm length of 25.1-27.8mm (Table 2).
The head is dorsoventrally very flattened.Pelage coloration is relatively variable, more or less golden red at the base of the dorsal fur, to dark brown near the tips of the dorsal hairs, and lighter golden brown on the underparts (Fig. 5B).The ears have a triangular shape, with broadly rounded tips.The tragus is short and blunt.The wing membranes are dark brown.The base of thumbs and soles of hind feet have fleshy pads (Fig. 5B).
The skull is small (GLS: 11.91-12.60mm), lightly built and very flat (Fig. 5B).The rostrum is short.The sagittal crest is absent.The lambdoid crests are well developed.The dental formula is I2/3 C1/1, P1/2, M3/3 = 32.The first upper incisor (I 2 ) is bicuspidate, with small cusps on cingulum.I 3 is unicuspidate, about half the height and crown area of I 2 .The upper canine has a posterior supplementary cusp.A diastema between I 3 and the upper canine is clearly visible.The protocones of M 1 and M 2 are welldeveloped.M 2 appreciably exceeds M 1 in width, and its width clearly exceeds its length.M 3 is relatively smaller and without a metastyle.The three lower incisors are tricuspidate.The first (PM 2 ) and second (PM 4 ) premolars are approximately equal in height and crown area (Fig. 5B).

Remarks
In northern Indochina, T. tonkinensis sp.nov.can be distinguished from T. fulvida and T. pygmaea by its significantly larger body and skull size (Table 2; Figs 4-5; see Feng et al. 2008 for comparisons with T. pygmaea), by K2P distances of at least 12% for Cytb and COI sequences and by K2P distances of at least 1.5% for the concatenation of the seven nuclear genes (5604 nt) (Appendix 5).Within the T. robustula complex, T. tonkinensis sp.nov. is morphologically overlapping with T. robustula, found in Sumatra, and T. malayana (= Tr2 haplogroup), collected from the Southeast Asian mainland, but differs from the first taxon by K2P distances of at least 5.2% in COI sequences (Appendix 5) and from the latter by K2P distances of at least 5.5%, 8.6% and 0.4% calculated from COI sequences (657 nt), Cytb sequences (1140 nt) and the concatenation of seven nuclear DNA sequences (5604 nt), respectively (Appendix 5).

Ecology and habitat
Like other species of Tylonycteris, T. tonkinensis sp.nov. is associated with woody bamboo groves.The new species is usually found in sympatry with the smaller species T. fulvida (Fig. 1).In northwestern Vietnam, bats of the new species were captured in mist-nests set near bamboo groves in forest edges adjacent to rural-residential areas at relatively high elevations, e.g., at 1010 m a.s.l in Hang Kia, Pa Co Nature Reserve (Hoa Binh Province) or at 1286 m a.s.l in its type locality.In Laos, the known localities are found at lower elevations, between 500 and 800 m a.s.l.

Distribution
Currently, the new species is known to occur in north-eastern Laos and northern Vietnam only (Fig. 1).

Cryptic species diversity in the genus Tylonycteris
Previous studies have detected high levels of genetic and karyological variation among specimens identified as either T. pachypus or T. robustula collected from different geographic locations in mainland Southeast Asia (Francis et al. 2010;Huang et al. 2014).Our COI analyses further revealed the existence of three divergent geographic haplogroups for both T. pachypus and T. robustula, for each corresponding to Sumatra, northern Indochina and southern Indochina (but also northwestern India and Peninsular Malaysia for T. robustula).Our Cytb dataset confirmed the existence of these geographic haplogroups in mainland Southeast Asia.In addition, a specimen of T. pachypus collected from Borneo and originally described as T. robustula was found to be highly divergent from the two other Indochinese haplogroups (Fig. 2).In the absence of Cytb data for Sumatran specimens, it is impossible to know whether they share the same mitochondrial lineage as those from Borneo.
Taken together, our mtDNA analyses show that all haplotypes sequenced for insular Tylonycteris are very divergent from those identified in mainland Southeast Asia.However, genetic inferences based on the maternally inherited mitochondrial genes are prone to be discordant with the true evolutionary history of the taxa, owing to various evolutionary processes, such as mtDNA introgression, incomplete lineage sorting or female philopatry (Avise 2000;Ballard & Whitlock 2004;Hassanin & Ropiquet 2007;Nesi et al. 2011;Rivers et al. 2005).Here, the geographic pattern of mtDNA diversity observed for the two species of Tylonycteris could be the consequence of female philopatry, i.e., the behavior of remaining in, or returning to the natal territory.Indeed, bat species with philopatric females generally display high geographic structure when relationships are examined with maternally inherited markers, such as the mitochondrial DNA.This pattern can disappear with biparentally inherited markers when adult males are able to disperse over long distances, allowing gene flow between otherwise isolated populations (Castella et al. 2001;Hassanin et al. 2015;Hulva et al. 2010;Pereira et al. 2009;Rivers et al. 2005).Behavioral and population genetic studies in southern China have shown that T. pachypus bats are philopatric to their natal area and that philopatry is especially pronounced in females (Hua et al. 2011(Hua et al. , 2013)).Although no data are available for T. robustula, female philopatry can also be predicted for this species, because it shares similar morphological, behavioural and ecological traits with T. pachypus (Medway 1972;Medway & Marshall 1970, 1972;Zhang et al. 2007).The social organization of T. pachypus and T. robustula, combined with their fragmented habitats, is therefore expected to result in limited gene flow between populations, especially among matrilines from distant geographic localities.For both species, this prediction is corroborated by the analyses of mtDNA markers, with the identification of three divergent geographically non-overlapping haplogroups.For T. robustula, this phylogeographic pattern is also supported by the nuclear sequence data, as the two Indochinese clades were recovered monophyletic with all the three introns containing enough nucleotide variation at the intra-specific level, i.e., CHPF2, HDAC1 and TUFM (Appendix 8).By contrast, our nuclear analyses (Fig. 3; Appendices 4, 8) do not support the reciprocal monophyly of the two Indochinese clades of T. pachypus, suggesting that gene flow was maintained by male dispersal or, alternatively, that their separation was too recent to be detected with our nuclear genes.
Nucleotide distances estimated from mtDNA genes between northern and southern Indochinese populations of T. robustula were more than twice those of T. pachypus (6.5% vs 2.8% in COI; 9.5% vs 2.8% in Cytb), indeed supporting a more recent divergence for the latter taxon, if we assume equal evolutionary rates.Similarly, the nuclear distances between the two Indochinese clades of T. robustula were between 0.41 and 0.56%, which is more than twice those calculated between Indochinese individuals of T. pachypus (0-0.2%;Appendix 5) and in the range of interspecific distances found in other groups of Laurasiatheria, such as fruit bats of the tribes Myonycterini (Nesi et al. 2013) and Scotonycterini (Hassanin et al. 2015), or cattle and buffalo of the tribe Bovini (Hassanin et al. 2013).Although none of the nuclear markers could be sequenced for Sumatran and Bornean Tylonycteris, their high mtDNA divergence from Indochinese populations (> 5.7 % in both COI and Cytb sequences; Appendix 5) suggests they might represent distinct lineages based on nuclear markers as well.In agreement with this view, our multivariate morphological analyses revealed that Indochinese bats of the T. pachypus complex constitute a distinct group separated from those collected on Sumatra.For the T. robustula complex, morphological overlap between haplogroups is more extensive, but pairwise comparisons of their PC mean scores support the distinctness of adjacent geographical taxa, such as Tr2 and Tr3 on the Southeast Asian mainland.
The close morphological similarity among taxa of Tylonycteris suggests that they have evolved under the influences of similar and specialized habitats, i.e., woody bamboo vegetation.Molecular evidence indicates, however, that T. pachypus should be split into at least two distinct species, T. pachypus on the Sunda islands (Sumatra and/or Borneo) and T. fulvida in mainland Southeast Asia, and that T. robustula should be divided into at least three species, with T. robustula on Sumatra, T. malayana in southern and western mainland Southeast Asia, and T. tonkinensis sp.nov. in northern Indochina.

The evolution of Tylonyteris spp. in Southeast Asia during the Pleistocene
Given that both species complexes, here named T. pachypus s. lat.and T. robustula s. lat., are usually found in sympatry across most of their geographic ranges in Southeast Asia, they are expected to share a common phylogeographical history.Our estimates of divergence times based on mtDNA sequences suggest that the genus Tylonycteris diversified during the Pliocene epoch (Cytb: 5.92 ± 0.65 Mya; COI: 4.56 ± 0.72 Mya) (Table 1).During the Miocene and until the early Pliocene, Southeast Asia was generally covered by large tracks of rain forests as a consequence of warm and humid climatic conditions (Meijaard & Groves 2006;Morley 2000).Thus, ancestors of both Tylonycteris species complexes were presumably widely distributed across Southeast Asia during the Pliocene.
Our molecular dating estimates indicate that the basal geographic splits within the two species complexes, i.e., between mainland Southeast Asia and Sumatra, took place approximately at the same time during the Early Pleistocene (between 2.70 and 1.96 Mya for T. pachypus, between 3.07 and 2.22 Mya for T. robustula; Table 1).The Pleistocene epoch is characterized by the onset of repeated cycles of cold glacial and warm interglacial periods as the results of the glaciations/deglaciations of the Northern Hemisphere, which implied contraction and expansion of rain forests in Asia (An et al. 2001;Meijaard & Groves 2006;Morley 2000).As bats of the genus Tylonycteris are highly dependent on woody bamboo vegetation for roosting, foraging and mating (Kunz 1982;Medway 1972;Medway & Marshall 1970, 1972), their Pleistocene biogeographic history was firmly constrained by the distribution of such bamboo habitats.The current distribution of woody bamboo species in Asia (Bystriakova et al. 2003b;Fig. 6) indicates that eight disjunctive biogeographic regions have higher species richness (> 5 species) for some bamboo genera: southern India (Ochlandra Thwaites), northern Myanmar (Cephalostachyum Munro), southern China (Dendrocalamus Nees and Bambusa Schreb.),Hainan Island (Bambusa), northwestern Thailand (Dendrocalamus and Gigantochloa Kurz ex Munro), Peninsular Malaysia (Gigantochloa), Sumatra and Borneo (Gigantochloa and Schizostachyum Nees).All these regions may therefore have acted as distinct bamboo forest refugia during the glacial periods of the Pleistocene (Fig. 6).Evidence for a number of these postulated glacial refugia has been reported in previous studies for many organisms, including bats (Flanders et al. 2011;Khan et al. 2010;Lin et al. 2014;Mao et al. 2013).Accordingly, we propose that the contraction of woody bamboo forests into different glacial refugia had fragmented the distribution of the Pliocene ancestors of both T. pachypus s. lat.and T. robustula s. lat.
In addition, we can assume that Pleistocene glacial periods resulted in higher interspecific competition between co-distributed species of Tylonycteris, because the supply of most suitable resources was more limited in glacial bamboo forest refugia (Medway & Marshall 1970).As noted in previous studies, T. pachypus s. lat.has a more manoeuvrable flight in cluttered habitats and forages on smaller insects than T. robustula s. lat.(Zhang et al. 2005(Zhang et al. , 2007)).Moreover, Medway & Marshall (1970) found that the smaller T. pachypus s. lat.can roost in the internodes with small entrance holes, which the larger T. robustula s. lat. is unable to enter.These differences suggest, therefore, that the smaller T. pachypus s. lat.have greater advantages than the larger T. robustula s. lat. in interspecific competition when natural resources are limited.Hence, isolated populations of T. robustula s. lat.may have been more exposed to bottlenecks and therefore more vulnerable to local extinction than those of co-distributed T. pachypus s. lat.
During interglacial periods of the Early Pleistocene, warmer and humid conditions resulted in the expansion of woody bamboo forests, which in turn may have favored the restoration of connectivity between isolated populations of both complexes.However, the isolated populations may have been connected or not, depending on their dispersal capacity and the distances between refugia.In T. robustula s. lat., these processes may have taken longer, because of its lower population abundance (Lande & Barrowclough 1987; Shaffer 1981), and may have been prevented in cases of extinction of transitional populations (Huang et al. 2014 and references therein).This scenario is supported by the fact that T. pachypus s. lat. is usually found to be more abundant than T. robustula s. lat. in bamboo forests (Zhang et al. 2004;Medway & Marshall 1972) and by the wider geographic range of T. pachypus s. lat.(Bates et al. 2008a(Bates et al. , 2008b;;Fig. 1).Knowing this, the body size differences between the two species complexes may be the key factor explaining why the basal divergence of northern Indochinese populations occurred earlier in T. robustula s. lat.(i.e., T. tonkinensis sp.nov) than in T. pachypus s. lat., i.e., 2.97-1.70vs 1.35-0.79Mya (Table 1).During Pleistocene interglacials, exchanges of Tylonycteris spp. between the continent and the islands of Sundaland were probably prevented because of the long distances between the glacial forest refugia, as well as the higher sea levels (Fig. 1).Our molecular dating estimates corroborate this scenario, as continental populations of Tylonycteris spp.from Indochina and Peninsular Malaysia diverged from insular populations (Sumatra and Borneo) in the Early Pleistocene (Table 1).

Implications for conservation
Previous studies considered T. pachypus and T. robustula to be common species and thus classified them as Least Concern in the IUCN Red List (Bates et al. 2008a(Bates et al. , 2008b)).Since our study reveals that both species in fact represent several species with more restricted distributions, the IUCN status of the different taxa should be reassessed urgently, including that of the new species, T. tonkinensis sp.nov.In addition, our study suggests that several biogeographic regions have acted as Pleistocene glacial refugia.This information is very important for developing more effective conservation strategies, particularly

Fig. 1 .
Fig. 1.A-B.Maps of Asia showing the distribution range (shaded) and type localities of described subspecies of T. pachypus (Temminck, 1840) and T. robustula Thomas, 1915 (Bates et al. 2008a, 2008b).C. Localities of Tylonycteris specimens included in this study.Triangles and circles refer to T. pachypus and T. robustula, respectively; the colours indicate the mtDNA haplogroups found in the Bayesian analyses of COI and Cytb sequences (see Fig. 2 for details).

Fig. 4 .
Fig. 4. Scatter plots obtained from morphological analyses of Tylonycteris spp. A. Range of GLS measurements of specimens within each group of Tylonycteris spp.B. Range of PC*1 scores of specimens of Tylonycteris spp.obtained from PCA of log-transformed raw data of craniodental measurements.C. Plot of PC 1 against PC 2 obtained from PCA on log-transformed standardized data.Triangles and circles refer to T. pachypus s. lat.and T. robustula s. lat., respectively.Colour patterns indicate the mtDNA haplogroups: green for Tp3 and Tr3 (northern Indochina); blue for Tp2 and Tr2 (other regions of the Southeast Asian mainland); red for Tp1 and Tr1 (Sundaland).

Fig. 6 .
Fig. 6.Distribution and species richness for the eight genera of woody bamboo forests currently found in Asia (adapted and modified from Bystriakova et al. 2003b) and putative geographic range areas of the species of Tylonycteris.During glacial periods of the Pleistocene, the highlighted biogeographic regions may have acted as bamboo refugia for Tylonycteris spp.Moreover, sea level falls (at around 120 m below present) exposed land bridges, such as the entire Sunda Shelf (adapted from Voris 2000), connecting the islands to the continent.

Table 2 .
Selected external and craniodental measurements (in mm) of Tylonycteris spp.Values are given as mean, standard deviation (SD), (n) and min-max.Acronyms and definitions for measurements are given in Material and methods.

Table 3 .
Factor loadings of craniodental characters for the PCs obtained from principle component analyses (PCAs) of specimens of Tylonycteris.Acronyms and definitions for measurements are given in Material and methods.The first PCA is based on raw data (PC*1), while the second PCA is based on the log-transformed, standardized data (PC1 and PC2).Values in bold indicate the most significant loadings.