Geographical Variation of Rhinolophus affinis (Chiroptera: Rhinolophidae) in the Sundaic Subregion of Southeast Asia, including the Malay Peninsula, Borneo and Sumatra

Rhinolophus affinis sensu lato is a widespread bat species in South and Southeast Asia which shows considerable geographical variation in its morphology, echolocation call frequencies and genetics. The taxonomic status of the taxon in the Sundaic subregion remains uncertain however as the limited studies to date have been largely based on morphology. The aim of the present study was to determine the taxonomic status of subspecific forms recognized in the subregion and to evaluate phylogeographic distinctiveness between those occurring in Borneo and the Malay Peninsula using genetic, morphological and acoustic datasets. Two forms were confirmed: R. a. nesites from Borneo and R. a. superans from the peninsula. The previous recognition of a population from southernmost Sumatra as R. a. superans was not supported, however, as this form is likely R. a. affinis. Genetic divergence between these three forms is rather deep and is estimated to have occurred during the arid climatic period of the Pleistocene when suitable habitats were reduced to isolated pockets. Our results support the phylogeographic distinctiveness hypothesis as R. affinis sensu lato shows discrete affinities between Borneo and the Malay Peninsula. Discovery of new forms of R. affinis is likely with greater sampling effort throughout the region. Our study also demonstrates the importance of employing multiple datasets in taxonomic evaluations, as the use of morphological and/or acoustic datasets alone could lead to erroneous conclusions.

The status of two subspecies, R. a. macrurus and R. a. superans, has recently been confirmed in continental Southeast Asia (Ith et al., 2015). The geographical boundary between these two forms lies in north Peninsular Thailand and accords with biogeographical demarcations within the region (Hughes et al., 2003(Hughes et al., , 2011de Bruyn et al., 2005;Woodruff and Tur ner, 2009). Rhinolophus a. macrurus, the Indo chinese form, exhibits considerable variation in its genetics, morphology and echolocation call frequencies (Ith et al., 2015). In contrast, the taxonomic status of the Sundaic form, R. a. superans, remains problematic, particularly in relation to populations on the island of Sumatra. Though Andersen (1905) described the Sumatra form as resembling specimens from the Malay Peninsula in cranial, dental and external morphology, the taxon has not been evaluated since this publication and its genetic and acoustic variation is unknown. Rhinolo phus a. superans is distributed throughout the Ma lay Peninsula (King sada et al., 2011;Ith et al., 2015), southern Sumatra (Huang et al., 2014) and central and north Sumatra (Andersen, 1907;van Strien, 1996;Csorba et al., 2003). The taxonomic status of R. a. nesites Ander sen has also not been evaluated. This form was proposed by Andersen (1905) as an offshoot of R. supe rans in Bunguran Island, north Natunas (ca. 230 km to the northwest of Borneo). The comparison was mainly based on the remaining parts of a damaged holotype which showed R. nesites has large ears, a broad horseshoe and a short tail. Though the form is recognized in recent literature (Med way, 1977;Koopman, 1994;Csorba et al., 2003;Sim mons, 2005), very little taxonomic work has been undertaken to confirm its status. The distribution of this subspecies from Bor neo includes Sabah, Sara wak and Kalimantan (e.g., Khan et al., 2008;Fran cis et al., 2010).
The use of multiple datasets strengthens the validity of taxonomic decisions. For instance, R. a. superans from the northernmost Malay Peninsula could be mistakenly assigned to R. a. macrurus on the basis of acoustic or morphological data alone, as this population has intermediate craniodental characters and a similar call frequency to R. a. ma cru rus but differs genetically (Ith et al., 2015). Similarly, the morphological cryptic Hipposideros bicolor (King ston et al., 2001) might not have been discovered without genetic and ecological data. Conver sely, genetics alone would not adequately discriminate the taxonomic status of other taxa such as R. ma crotis and R. siamensis as these show very shallow genetic differences (Francis et al., 2010). Simi lar cases include Miniopterus schreibersii (Fur man et al., 2010), Eptesicus serotinus, E. nilssonii (Mayer and von Helversen, 2001) and Myotis annamiticus (Kruskop and Tsytsulina, 2001;Francis et al., 2010).
Rhinolophus a. superans may have similar morphological and genetic variation to that found in the Indochinese form of R. affinis: R. a. macrurus (Ith et al., 2015). Francis et al. (2010) have shown that widespread taxa often have substantial geographic variation in their barcode sequences and that populations from Peninsular Malaysia and Borneo are often genetically distinct (e.g., Khan et al., 2008Khan et al., , 2010. As a consequence, the aim of the current study was to determine the taxonomic status of R. a. superans and R. a. nesites and evaluate the phy lo geographic distinctiveness of R. affinis from Bor neo and the Malay Peninsula using a combination of genetic, morphological and acoustic datasets.

Study Specimens and Sampling Sites
Seventy-six specimens were available for morphological study, including five from south-western Sumatra, seven from Sarawak, north-western Borneo and 64 from the Malay Peninsula. Two specimens from Central Java, Indonesia and two specimens from Musoorie, northern India were also included for comparison. Samples examined were from existing museum collections and those arising from recent surveys. Specimens were examined in collections held at the Princess Maha Chakri Sirindhorn Natural History Museum, Prince of Songkla Uni versity, Thailand (PSU collection); Harrison Institute, UK (HZM collection); Museum Zoologicum Bogoriense, Research Center for Biology-Indonesian Institute of Science, Indonesia (MZB collection); Museum of Texas Tech University, USA (TTU collection); and Zoological Museum of Universiti Malaysia Sarawak, Malaysia (UNIMAS collection).
Specimens from the Malay Peninsula were collected by Saveng Ith and the Small Mammals and Birds Research Unit Team of PSU between August 2011 and May 2012. Bats were surveyed in the field using a combination of harp traps, mist nets and hand nets and were captured and handled in accordance with guidelines approved by the American Society of Mam malogists (Gannon et al., 2007). Field surveys were conducted in several localities in Thailand including Hala Bala Wildlife Research Station, Khao Namkhang National Park, Khao Ban Tad Wildlife Sanctuary, Rajjaprabha Dam and Ton Nga Chang Wildlife Sanctuary. All study localities where the 76 specimens were collected are illustrated in Fig. 1

Morphological Measurements
Thirty-three external and craniodental characters of each specimen were measured following Bates and Harrison (1997), Thomas (1997), Csorba et al. (2003) and Furey et al. (2009). External characters were measured using a pair of dial calipers to the nearest 0.1 mm and craniodental characters were measured to the nearest 0.01 mm using a digital caliper under a stereo microscope. Definitions for external measurements are as follows, FA: forearm length -from the extremity of the elbow to the extremity of the carpus with the wings folded; EL: ear length -from the lower border of the external auditory meatus to the tip of the pinna; TL: tail length -from the tip of the tail to its base adjacent to the anus; HF: from the extremity of the heel behind the os calcis to the extremity of the longest digit, not including the hairs or claws; TIB: tibia length -from the knee joint to the extremity of the heel behind the os calcis; 2MT, 3MT, 4MT, 5MT: length of metacarpals -taken from the extremity of the carpus to the distal extremity of the second, third, fourth and fifth metacarpals, respectively; 1P3D, 2P3D, 1P4D, 2P4D, 1P5D, 2P5D -length of the first and second phalanges of the third, fourth and fifth digits, respectively -taken from the proximal to the distal end of the phalanx; GWN -greatest width of noseleaf -greatest diameter across the horseshoe; GHN: greatest height of noseleaf -from the base of the horseshoe to the tip of the lancet, not including the hairs.
All skulls were extracted for examination. Definitions for craniodental measurements were as follows: SL: skull lengththe greatest length from the occiput to the front of the canine; CCL: condylo-canine length -from the exoccipital condyle to the anterior alveolus of the canine; ALSW: the greatest width across the anterior lateral compartments of the rostrum; AMSW: anterior median swellings width -the greatest width across the median swellings in dorsal view; ZYW: zygomatic width -the greatest width of the skull across the zygomata; BW: braincase width -width of the braincase at the posterior roots of the zygomatic arches; GBW: greatest braincase width -width of the braincase, the greatest width across the braincase; BOW: basioccipital width -least distance between the cochleae; MAW: mastoid width -greatest width of the braincase taken across the mastoid region; IOW: interorbital width -the narrowest width of interorbital constriction; PB: palatal bridge -length of bony palate excluding the posterior spike; M 3 M 3 W: posterior palatal width -taken across the widest part of the outer borders of the third upper molar; C 1 C 1 W: anterior palatal length - Dashed arrows indicate type localities and subspecies names, solid arrows indicate the transition zone of biota within the Malay Peninsula, dashed lines indicate the echolocation frequencies (min-max), and the two-headed arrows indicate the echolocation frequencies (min-max) as a whole from each echolocation zone. Note: the northern boundary of the Sundaic subregion is sometimes placed at the Isthmus of Kra (e.g., Lekagul andMcNeely, 1988 andCorbet andHill, 1992) taken across the widest part of the outer border of the upper canine; CM 3 L: upper toothrow length -from the front of the upper canine to the back of the crown of the third upper molar; CM 3 L: lower toothrow length -from the front of the lower canine to the back of the crown of the third lower molar; ML: mandible length -from the most posterior part of the condyle to the most anterior part of the mandible, including the lower incisors; CPH: least height of the coronoid process -from the tip of the coronoid process to the apex of the indentation on the inferior surface of the ramus adjacent to the angular process. Baculum characters were measured to the nearest 0.01 mm using a digital caliper under a stereo microscope. Thirty bacula were available for examination, comprising 27 from the Malay Peninsula, two from Sumatra and one from Borneo.

Echolocation Call Measurements
Values for the frequency of maximum energy (FMAXE) for R. affinis in this study were obtained from field work. To avoid pseudo-replication, one echolocation call per bat was used in analysis. In total, 72 calls (from 72 bats) were available for meas urement. Fifty-nine calls were from the Malay Peninsula, one from north-western Borneo, six from Central Java and five from southwestern Sumatra.
Echolocation calls were recorded from bats held in the hand using a Pettersson D-240X bat detector and in some instances, a Pettersson D1000X bat detector (Pettersson Elektronik, AB). The Pettersson D-240X detector was set in ×10 time-expansion mode and call data was recorded to a digital iRiver iHP-120 Multi Codec Jukebox recorder. Where a Pettersson D1000X was used, calls were stored on a built in Compact Flash (CF) card (type I). The detector was set to manual recording mode (MAN) and the maximum sampling frequency (fs) to 768 kHz. A time expansion factor of ×10 was also used. All sound files were recorded and saved in 'wav' format for analysis. Call components were displayed using spectrogram, oscillograms and power spectrums in BatSound Pro 3.31 (Pettersson Elektronik, AB) in which sampling frequency was formatted as 44.10 kHz and spec trograms were set to 1,024 sampling size using Fast Fourier Transforms (FFT) with Hanning windows. In all cases, FMAXE (kHz) was measured from the constant frequency portion of a call using power spectra and the mean value was used in analysis.

Morphological and Acoustic Analyses
Statistical analyses were carried out using SPSS 16.0 (SPSS Inc., Chicago, USA) and PC-ORD 5.10 for Windows (MjM Software, Oregon, USA). Descriptive statistics (minimum, maximum, mean and standard deviation) were calculated for echolocation, external and craniodental measurements. Normal ity of data and homogeneity of variances were tested prior to using parametric t-tests and non-parametric Mann-Whitney U-tests to evaluate sexual dimorphism in size. Multiple comparisons of char acters between populations were calculated us ing a multivariate analysis of variance (MANOVA). Principal component analysis (PCA) on the correlation matrix was used for multivariate comparisons.

Molecular Analysis
Tissue was collected from different organs of voucher specimens such as liver, tongue and wing membrane and preserved cold in 95% concentration ethanol. Two mitochondrial DNA gene fragments were used for phylogenetic analysis. A 657 base pairs segment of 17 sequences of cyto chrome c oxidase I (COI) was analyzed at the Canadian Center for DNA Barcoding (CCDB) using the barcoding protocols, methods of analyses were detailed in Ivanova et al. (2012). A 832 base pairs segment of 19 sequences of cytochrome b (Cytb) gene was generated and analyzed in collaboration with the Coral Triangle Partnerships in Inter national Research and Educa tion Project (https://sci.odu.edu/ impa/ctpire.html). Ge nom ic DNA was isolated from bat tissue samples using the Qiagen DNeasy mini kit (Qiagen, Valencia, CA) following manufacturer's instructions and Cytb sequences were generated, aligned and proofread as described in Willette and Padin (2014) using the primers Cytb 07 (5'-AATAGGAGG-TATCATTCGGGT-3') and Cytb 09 (5'-GTGACTTGAAAA-ACCACCGTT-3'). The full lengths (1,140 base pairs) of 13 Cytb sequences and 413 base pairs segment of five Cytb sequences were analyzed (DNA extraction, PCR amplifications, and sequencing reaction) by F.A.A.K. following Khan et al. (2013) using primer set LGL765 (5'-GAAAAACCAYCGTTG-TWATTCAACT-3'), LGL766 (5'-GTTTAATTAGAATYTYAG CTTTGGG-3') with an annealing temperature of 50˚C.
In total, 37 Cytb sequences and 17 sequences of COI were available. Sequences from GenBank and Barcode of Life Data Systems (BOLD) were also accessed, and eight sequences of Cytb gene (accession number: EF108156-EF108160, EU521607, JN106274 and JN106280) from Borneo and Peninsular Malaysia were included for comparison. Twenty-one sequen ces of COI gene were included, 11 se quences were from Peninsular Malaysia (accession no: HM541330-HM541332, HM541407-HM541414) and 10 sequences from Peninsular Thailand.
Phylogenetic relationships among sequences were reconstructed using maximum-likelihood in the MEGA 5.2.2 program (Tamura et al., 2011). The most appropriate substitution model was determined using Akaike Information Criterion (AIC) and Bayesian Information Criterion (Bickham et al., 2004) as implemented in jModelTest 2.14 (Darriba et al., 2012). Among the 88 models in the 100% confidence interval, the Hasegawa-Kishino-Yano substitution model (HKY) was the best-fit model selected for COI. While General Time Reversible models (GTR) with gamma distribution (G) were the best-fit model selected for Cytb. We also performed Bayesian Analysis using MrBayes 3.2.2 (Huelsenbeck and Ronquist, 2001). In Bayesian Analysis, convergence stationary was search ed by two independent Metropolis-coupled Markov chain Monte Carlo (MCMC), each comprising three incrementally heated chains and one cold chain, run for 24 million generations, with parameters sampled every 500 generations. Con vergence stationary of the MCMC chains was evaluated by inspecting whether the standard deviation of split frequencies < 0.01 and the potential scale reduction factor (PSRF) reached 1.0 for all parameters. We also investigated the convergence using Tracer 1.5 . 12,000 trees of initial phase of the Markov chain were discarded as 25% burn-in. A congeneric Rhi nolophus stheno was used as an out-group in the phylogenetic analysis of Cytb gene in order to examine the monophyletic lineage of R. affinis.
To estimate the time to the most recent common ancestor (TMRCA) among the observed clades, Cytb gene was analyzed in BEAST 1.8 . GTR + G was selected as the best substitution model based on jModelTest and a relaxed-clock model with an uncorrelated lognormal distribution was selected for the substitution rate. We performed two independent runs of MCMC chains with 60 million generations with parameters logged every 1,000 generations. Tracer 1.5  was used to combine the two runs as well as to examine the effective sample size (ESS) for the parameters. Trees were collated using Tree Annotator 1.8 where Maximum clade credibility tree and Median heights were selected; and 10% (6,000 trees) sample trees were selected as burn in. To convert the estimates scaled by mutation rate to calendar years, we used the mean substitution rate of 1.30 × 10-8 subs/site/year which was previously used in hipposiderid bats (Thong et al., 2012;Lin et al., 2014). To calculate the genetic distance within and between clades, pairwise genetic distances (P-distance model) in MEGA 5.2.2 were computed.

Morphology
To explore sexual dimorphism, localities where male and female specimens were both available were selected; numbers of each were adjusted to balance sample sizes and so 22 males and 22 females were compared. No significant differences were found in 33 external and cranial characters between the sexes. A total of 12 external and cranial characters were retained for multivariate analysis, these being selected on the basis of their eigenvector values in a preliminary PCA. A PCA using these 12 characters for 74 specimens from the Sundaic subregion generated four relatively isolated groups including Borneo, Sumatra, southern Malay Peninsula and northern Malay Peninsula groups (Fig. 2). Spec imens from Borneo exhibited a higher degree of isolation among the groups.
Specimens from Borneo were distinguished from Sumatra and Malay Peninsula specimens by their generally smaller external and cranial measurements and noseleaves. Specifically, Borneo specimens were smaller on average in FA, TL, TIB and HF (P < 0.05) and several wing measurements (2MT, 3MT, 4MT, 5MT and 1PD3; all P < 0.05). Several skull characters were also significantly smaller, including SL, ZYW, CM 3 L, C 1 C 1 W, M 3 M 3 W, CM 3 L and CPH (all P < 0.05 - Table 1). The skull of these specimens has a short frontal depression and the canines and other teeth are smallest overall (Fig. 3).
The noseleaf is small, as is GWN with an average width of 9.1 mm, while GHN is also small, at 12.9 mm. The median emargination of the horseshoe is narrow (Fig. 4C). The rudimentary secondary noseleaf is less developed and completely concealed by the horseshoe and surrounding dense hair (Fig.  4C). The sella is small and slender, rounded off on the top and the lateral margin is more strongly constricted in the middle (Fig. 5C). The internarial cup is moderate in size and the margin is developed (Fig.  4C). The connecting process is small, slender, rather pointed and covered with numerous short hairs and shows the notch pattern on the top. The lancet is small, slender, triangular-shaped and straight-sided.
Specimens from Sumatra also formed a relatively isolated group (Fig. 2). Compared with specimens from the northern Malay Peninsula, Sumatran specimens are externally smaller in TIB, 2P3D, 1P4D and 2P5D (P < 0.05) but larger in  4D and Table 1). The skulls of Sumatran specimens also have significantly smaller MAW, GBW, ALSW, AMSW, IOW, CM 3 L and CM 3 L (all P < 0.05 -Figs. 3, 6, 7, and Table 1). Compared with specimens from the southern Malay Peninsula, Sumatran specimens are similar in size with only two external (TIB and P2D5) and one craniodental character (AMSW) significantly smaller, and three characters significantly larger (CCL, PB, and C 1 C 1 W -P < 0.05). Sumatran specimens were found to have a more developed sagittal crest (Fig. 3B) however, which is well built and visible from the supraorbital ridges to the lambda. The noseleaf of Sumatran specimens is medium sized in general and shares many characters with specimens from the Malay Peninsula. GWN in Sumatran specimens is slightly smaller than Malay Peninsula specimens with an average of 9.9 mm; GHN is highest in the Sumatran population, averaging 15.0 mm. The median emargination of the horseshoe is as wide as specimens from Central Java and Malay Peninsula and differs from specimens from Sarawak and India (Fig. 4D). The rudimentary secondary noseleaf is visible in dorsal view, with fewer hairs compared to Sarawak and Central Java specimens (Fig. 4D). The sella is large, tall and rounded off on the top, and the lateral margin is only slightly constricted in the middle (Fig. 5D). The internarial cup is moderate in size and the margin is less developed compared to specimens from Sara wak (Fig. 4D versus Fig. 4C). The connecting process is typically round and the lancet is triangular, straight-sided and high.
Specimens from the Malay Peninsula had the largest craniodental measurements overall (Table 1). The rostral chambers are large (Fig. 6D) and ALSW and AMSW are broad, averaging 6.15 mm and 4.26 mm, respectively. The anterior median swellings are inflated (Fig. 3D) and rounded in the dorsal view (Fig. 6D). The frontal depression (Fig. 3D) and supraorbital ridges (Fig. 6D) are elongated and the palatal bridge is long, with CM 3 L, ML and CM 3 L also large (Fig. 7D). Similarly, the noseleaf is relatively large with the largest GWN, averaging 10.0 mm. The rudimentary secondary noseleaf is developed but almost invisible in the dorsal view being largely concealed by the horseshoe (Fig. 4E,  4F). The sella is very broad and lacks an obvious middle constriction as the lateral margins gradually constrict (Fig. 5E, 5F). The tip of the sella is always rounded off. The internarial cup is broad with welldefined but not especially developed lateral margins 3.16 ± 0.14 ( Fig. 4E, 4F). The connecting process is typical of the species, large and rounded off and covered with many short hairs. The lancet is broad and high with elongate tip, and its lateral margins are normally straight-sided or slightly convex at the base in some individuals.

Baculum
The bacula of Sumatran specimens (n = 2) are similar to those from the Malay Peninsula, although some differences are apparent. Overall, the bacula of Sumatran specimens are slightly shorter and the basal portion is more inflated and rounded (Fig. 8B  versus 8C). In the lateral view, the bacula of Su matran specimens have a larger shaft and an enlarged and less pointed tip. An enlarged tip is also found in many but not all Malay Peninsula specimens. In the dorsal view, the vertical ridges along either side of the basal part are almost invisible and sometimes absent in Sumatra specimens but are well developed in Malay Peninsula specimens.
The baculum of the specimen from Sarawak is similar to those of Sumatran specimens, just slightly more slender overall with a less inflated basal portion ( Fig. 8A versus Fig. 8B). In the lateral view, the tip portion of the shaft is also swollen in character but is elongated and less prominent (Fig. 8A) compared to Sumatran specimens (Fig. 8B). In the dorsal view, the basal emargination is deeper and narrower.

Echolocation
Extensive variation in call frequencies occur within the Sundaic subregion, with differences of 20 kHz recorded across the range (62.3-82.3 kHz). Average frequencies observed are: Central Java

Genetics
Results from both maximum likelihood (ML) and Bayesian analysis (BA) showed similar topologies in phylogenetic trees. Three clades were recovered based on Cytb genes (Fig. 9). Clade A and C lineages were supported by high bootstrap values (BT = 90-99%) and posterior probabilities (PP = 100%) while clade B was supported by lower BT = 60% but rather high PP = 94%. The recovery of the three lineages was very consistent in the analyses; however the recovery of basal linage was inconsistent. The two possible basal lineage relationships through our analyses (A and B, or B and C - Fig. 9) were poorly supported (e.g., BT = 30%, PP = 75%).
Clade A comprised sequences from Borneo, whereas clade B comprised sequences from Borneo, Central Java and Sumatra, and clade C comprised sequences from the Malay Peninsula (Fig. 10). Pairwise genetic distances within clades were low at 0.01%, 0.00-0.03 (mean, range) for clade A, 0.06%, 0.00-1.30 for clade B and 0.05%, 0.00-0.10 for clade C. Mean genetic distance between Borneo and Central Java-Sumatra was lower (clade A versus B: 2.8%, 2.6-3.3), and relatively higher between the Malay Peninsula and Borneo (clade C versus A: 3.7%, 3.7-4.4) and the Malay Peninsula and Central Java (clade C versus B: 3.6%, 3.0-4.4). Based on the mean genetic distance, the Central Java and Borneo clades (B and A) shared a more recent ancestor than the Malay Peninsula clade (C). Clade C was therefore assumed to be basal to clade A and B. Results from both ML and BA illustrated similar topologies, with two clades recovered for COI gene (Fig. 11). Clade A (BT = 99%, PP = 100%) comprised all sequences from the Malay Peninsula whilst clade B (BT = 59%, PP = 100%) comprised sequences from Sumatra (Fig. 12). Pairwise genetic distance within clades were low at 0.02%, 0.00-0.07 (mean, range) for clade A and 0.03%, 0.00-0.05 for clade B. Mean genetic distance between the clades was 2.2%, 1.7-2.7 (A versus B). As both clades were consistently recovered with strongly supported values and observed in Cytb analysis (clades B and C - Fig. 9), these populations were recognized as two isolated lineages.
Bayesian estimates of time to the most recent common ancestor (TMRCA) provided effective sample size (ESS) values > 500 for all parameters. The inferred TMRCA for all recovered clades, including the Borneo and Central Java and Malay Peninsula clades (A versus B, C) was 1.7 million years before present (Myr BP) (95% CI 1.09-2.35) (Fig. 9), corresponding to an early stage of the Pleisto cene epoch. The TMRCA for B versus C was more recent at 1.30 Myr BP (95% CI 0.82-1.81) However, as recovery of basal lineages was inconsistent (switching between clade A and C), we assume TMRCA between lineages is more or less the same (ca. 1.30-1.70 Myr BP).

Variation within the Malay Peninsula
Intraspecific variation was also found within the Malay Peninsula. Specimens from the high call frequency zone (green shading A: 77.3-79.3 kHz;   7D-E). However, both populations have similar bacula morphology. A PCA using nine external and cranial characters of all specimens from Malay Peninsula illustrated two relatively isolated groups (Fig. 13).

DISCUSSION
On the basis of morphology, bacula and genetic evidence, three geographical forms of R. affinis are recognized in the Sundaic subregion of South east Asia. Two of these are referred to their existing names (R. a. nesites from Borneo and R. a. superans from the Malay Peninsula), while the population from Sumatra is provisionally referred to R. cf. affinis due to its morphological and genetic differences from R. a. superans in the Malay Penin sula. Although sampling sizes for this regionally widespread species are limited, each genetic clade identified here corresponds to a unique morphology that reflects the phylogeographic distinctiveness of different locations. Similar divergence patterns have been found in other bat species from Borneo and the peninsula (Francis et al., 2010;Khan et al., 2010), and also in murine rodents (Gorog et al., 2004). Call frequencies in the region are not congruent with this pattern however and disparities between acoustic and other datasets have also been observed in R. af fi nis from the Indochinese subregion (Ith et al., 2015), as well as R. malayanus (Soi sook et al., 2008) and Hipposideros larvatus (Tha bah et al., 2006). Rhinolophus a. nesites was described as occurring from Sarawak (Bau, Kuap) to West and South Borneo (Medway, 1977) and this form was recognized by Koopman (1994) and Csorba et al. (2003). The holotype is deposited in the American Museum of Natural History (AMNH.104753♀) and is badly damaged, with only the teeth and lower jaw remaining in good condition. Based on this, Andersen (1905) described R. a. nesites as comparable to R. a. superans, but with a shorter TIB and smaller MT. Specimens subsequently collected from Sarawak agree with the holotype description in having a short TIB and MT (Table 1). However, they differ in having relatively smaller EL and GWN. This is probably because the holotype was described from Bunguran Island, and our data indicate that specimens from islands (e.g., Adang, Rawee and Koh  Fig. 1) were smaller in many instances compared to spec imens from northwards of Khao Namkhang National Park (T15) (the lower call frequency zone, orange shading B: 69.5-72.6 kHz - Fig. 1) particularly in cranial characters (Table 1). The former Surin) tend to have larger ears and noseleaf characters (e.g., horseshoe, sella and connecting process) and emit lower call frequencies (Table 2). Our comparisons also show that R. a. nesites is significantly smaller on average (P < 0.05, table of comparisons not included) than R. a. superans and R. cf. affinis (from southwestern Sumatra) in other external characters and many cranial characters. Rhinolophus a. nesites also differs in noseleaf and baculum characters and genetic data support this divergence as none of the sequences from Borneo (clade A) nested with Malay Peninsula sequences (clade C) or vice versa (Figs. 9 and 10). Andersen (1905) included Sumatra in the distribution of R. a. superans based on a specimen from Sirambas, central Sumatra, the only specific locality record from the region (Andersen, 1907). This was accepted by Csorba et al. (2003). In our study, R. cf. affinis from southwestern Sumatra differed in many skull and baculum characters from Malay Peninsula specimens and also genetically (COI and Cytb -Figs. 9 and 11). Sumatran specimens are more similar to Central Javan specimens craniodentally (Table 1) and genetically (Fig. 9). We therefore distinguish the southwestern Sumatran population from peninsular populations. However, since our sample was limited to the southern tip of Sumatra, the possibility that specimens from central and northern parts of the island could be allied with peninsular populations cannot be excluded as morphological and genetic affinities between adjacent areas of different islands have been found for R. affinis in Wallacea (Mahara da tunkamsi et al., 2000). The presence of two Cytb sequences (GenBank) in clade B (Fig. 9) from an unspecified locality or localities in Borneo requires comment. If these FIG. 9. Baysian phylogenetic tree based on Cytb gene. Scores on the branches refer to bootstrap support values (1,000 iterations) derived from maximum likelihood (1st score) and Bayesian posterior probabilities (2nd score); --= no support value. The horizontal bars on the tree branches represent the 95% highest posterior density intervals for the divergence estimates. Specimens are labeled by specimen codes (CHGTK, EF, EU, JN, IS, MZB, PS and TK) and collecting localities. The symbols of clades correspond to the genetic distribution map, Fig. 10 Millions of years FIG. 10. Distribution of Cytb clades of R. affinis within the Sundaic subregion. The shape of the symbols corresponds to clades defined in Fig. 9. Black symbols are sequences from the current study whereas grey symbols are sequences from GenBank. Localities of sequences not listed in the methods and materials of the current study are listed for the first time as following, BA = Jambusan Cave, Bau, Sarawak; GB = Gunung Berumput, Sarawak; GG = Gunung Gading NP, Sarawak; KM = Kabumen, Central Java and PC = Prachuap Kiri Khan. unK = unknown specific locality from Borneo (sequences from GenBank). Dashed arrows indicate the type localities of subspecies. Black solid arrows indicate the transition zones of biota in the peninsula sequences are genuinely derived from Borneo, then their true localities are most like ly to be from regions near to Java or Sumatra and result from dispersal during the Pleistocene. Another possibility is that R. a. affinis and R. a. nesites are sympatric in Borneo which would challenge their subspecieslevel status. Alternatively, the location of Borneo given for these sequences could be due to incorrect labeling, and they are from Java or Sumatra but represent a genetically distinct group. Rhinolophus a. superans was described as similar to the Indochinese form R. a. macrurus, but with a shorter TL, ALSW, AMSW, CM 3 L, CM 3 L and BW (Andersen, 1905). Ith et al. (2015) found that R. a. su perans is significantly smaller in many external measurements, larger in many skull measurements and distinguishable genetically (D-loop, COI). Our study indicates that morphological variation of R. a. superans in the Malay Peninsula aligns closely with the Kangar-Pattani Line (van Steenis, 1950;Whit more, 1984) and climatic zones (Hughes et al., 2011), although this was not supported genetically. Morphological differences were also observed between specimens from the mainland and adjacent islands (Adang, Rawee and Kho Surin) and as noted above, island specimens have larger noseleaves, rostral chambers and other skull measurements and also emit lower call frequencies.
These differences reflect known relationships between call frequencies and horseshoe bat skull characters (Heller and von Helversen, 1989;Francis and Habersetzer, 1998;Barc lay et al., 1999;Guillén et al., 2000), but were not supported genetically as all sequences of each gene nested together in one clade (clade C, Fig. 9 and clade I, Fig. 11).
The extensive geographical variation in echolocation call frequencies emitted by R. affinis sensu lato does not accord with morphological and genetic variation. Similar variation in call frequencies has also observed in R. a. macrurus (Ith et al., 2015) and R. malayanus (Soisook et al., 2008). In the Malay Peninsula, call frequencies for R. affinis from islands (loc. T4, T17 and M2) are lower than adjacent localities on the mainland. Two clear frequency patterns occur that correspond with morphological variation (Fig. 1). Higher frequency calls occur FIG. 11. Maximum likelihood tree based on COI. Scores on the branches refer to bootstrap support values (1,000 iterations) derived from maximum likelihood (1st score) and Bayesian posterior probabilities (2nd score); --= no support value. Specimens are labeled by specimen codes (IS, HM, HZM and MZB) and collecting localities. The symbols of clades correspond to the genetic distribution map, Fig. 12 FIG. 12. Distribution of COI clades of R. affinis in the Sundaic subregion. The shape of the symbols corresponds to clades defined in Fig. 11. Black symbols are sequences from the current study and Ith et al. (in review) whereas grey symbols are sequences from GenBank. Localities of the sequences not listed in the methods and materials of the current study are listed for the first time as following, ER = Endau Rompin National Park, Peninsular Malaysia; KL = Kuala Lompat, Pahang; NS = Negeri-Sembilan; TT = Thaninthary Div, Myanmar. Dashed arrows indicate the type localities of subspecies. Black solid arrows indicate the transition zones of biota in the peninsula south of the Kangar-Pattani Line in tropical evergreen rain forest and lower frequency calls occur north of the line in semi-evergreen rain forest. These differences are reflected in some external and craniodental characters (Fig. 13), but not by Cytb and COI genes ( Figs. 9 and 11). This suggests that these differences are not a result of selection for clear social communication (which would lead to phylogenetic distinctiveness: Heller and von Helversen, 1989;Kingston et al., 2000;Kingston and Rossiter, 2004), but may have recently evolved in association with climatic conditions, foraging habitats and/or prey availability. The current forest and climatic conditions of the region began as recently as the end of the last glacial maxima (ca. 9,500 years ago, after the breakup of Sunda Shelf land-bridge - Voris, 2000;Inger and Voris, 2001), which may not have been long enough for significant genetic differences to evolve. Our results suggest that R. a. nesites, R. a. superans and R. cf. affinis diverged in the early Plei stocene epoch (1.7-1.3 Myr BP). This may have been caused by refugial isolation prior to the coldest Pleistocene period and accords with estimated divergence times for other taxa in the region e.g. gymnures (Ruedi and Fumagalli, 1996), murine rodents (Gorog et al., 2004), bats (Khan et al., 2010;Lin et al., 2014), herpetofauna (Inger and Voris, 2001) and termites (Gathorne-Hardy et al., 2002). The formation of R. affinis lineages on the Sunda shelf may be partly explained by its ecology. Rhinolophus affinis is a cave-dwelling bat species which forages in the understorey of forest, including mature lowland rainforest, dry forest and disturbed areas (Francis, 2008). As such, the historical transition from a relatively stable tropical environment and perhumid climate during the Miocene (Gorog et al., 2004) to more arid and cool climatic conditions in the Plio-Pleistocene when suitable habitats in Southeast Asia were reduced to isolated pockets (Heaney, 1991;Morley, 1998Morley, , 2000van der Kaars et al., 2001) may explain the current biogeography of R. affinis.
In conclusion, our study demonstrates the importance of employing multiple datasets in taxonomic evaluations, as use of morphological and/or acoustic datasets alone could lead to erroneous conclusions. The discovery of additional population structures (e.g., R. cf. affinis from Sumatra) is also predicted in Southeast Asia with greater sampling effort throughout the region.