Mitochondrial genetic differentiation and morphological difference of Miniopterus fuliginosus and Miniopterus magnater in China and Vietnam

Because of its complicated systematics, the bent-winged bat is one of the most frequently studied bat species groups. In China, two morphologically similar bent-winged bat species, Miniopterus fuliginosus and Miniopterus magnater were identified, but their distribution range and genetic differentiation are largely unexplored. In this study, we applied DNA bar codes and two other mitochondrial DNA genes including morphological parameters to determine the phylogeny, genetic differentiation, spatial distribution, and morphological difference of the M. fuliginosus and M. magnater sampled from China and one site in Vietnam. Mitochondrial DNA gene genealogies revealed two monophyletic lineages throughout the Tropic of Cancer. According to DNA bar code divergences, one is M. fuliginosus corresponding to the Chinese mainland and the other is M. magnater corresponding to tropical regions including Hainan and Guangdong provinces of China and Vietnam. Their most recent common ancestor was dated to the early stage of the Quaternary glacial period (ca. 2.26 million years ago [Ma] on the basis of D-loop data, and ca. 1.69–2.37 Ma according to ND2). A population expansion event was inferred for populations of M. fuliginosus at 0.14 Ma. The two species probably arose in separate Pleistocene refugia under different climate zones. They significantly differed in forearm length, maxillary third molar width, and greatest length of the skull.


Introduction
Bent-winged bats Miniopterus schreibersii were previously considered to be widely distributed across the Old World ranging from Europe through to the Pacific (Koopman 1994;Wilson and Reeder 2005). The complex is now known to comprise several species (Appleton et al. 2004;Tian et al. 2004;Furman et al. 2009Furman et al. , 2010. In China, two bent-winged bat species, Miniopterus fuliginosus Hodgson, 1835 and Miniopterus magnater Sanborn, 1931 ( Fig. 1), were formerly confused with M. schreibersii (Hendrichsen et al. 2001;Appleton et al. 2004;Tian et al. 2004). Tian et al. (2004) argued that M. schreibersii from Guangxi and Hainan in China should be considered as M. fuliginosus based on intraspecific mtDNA divergence levels (Tian et al. 2004), which supported Maeda's recognition of Asian Miniopterus schreibersii as a distinct species, M. fuliginosus (Maeda 1982). Furthermore, Maeda regarded the bent-winged bats from Hainan Island to differ from individuals on the Chinese mainland (Maeda 1982).
Those specimens originally labeled as M. schreibersii in South-East Asia may represent M. magnater (Hendrichsen et al. 2001). Miniopterus magnater has been recorded in southern China, including Hainan, Guangdong, Hongkong, and Fujian (Smith and Xie 2008) and shares similar morphological characteristics with M. fuliginosus in China (Maeda 1982;Smith and Xie 2008). Most of their morphological characters overlap, and the two species differ only in skull size, with M. magnater being slightly larger and wider than that of M. fuliginosus. Confusion of morphological characteristics between M. fuliginosus and M. magnater makes the mapping of the distribution extremely difficult in China, which suggest that bent-winged bat species distribution limits and characters need further research and evaluation.
Molecular data could play a major role in a re-examination of the taxonomics, phylogeny, and lineage divergences of the bent-winged bats. More recently, cytochrome c oxidase subunit I (COI) has been selected as the DNA bar codes for members of the animal kingdom (Hebert et al. 2003;Marshall 2005;Kerr et al. 2009) and has established a standardized approach to help field researchers in identifying species accurate (Borisenko et al. 2008). And COI bar codes are proved as an effective tool for both differentiating and identifying species of bats (Borisenko et al. 2008;Francis et al. 2010). Francis et al. (2010) have sequenced the COI bar codes of M. fuliginosus and M. magnater (Francis et al. 2010), which are available for clarifying the taxonomy and distribution range of the two species in China.
Glaciation and resultant geographic isolation might be considered as major mechanisms underlying the genetic differentiation of these two closely related bat species, M. fuliginosus and M. magnater. As a particular region, China encompasses both the Palearctic and Oriental biogeographic regions and occupies several climatic zones. Due to the influence of Pleistocene glacial cycles, many mammals, including bats, show distinct genetic differentiation among populations from different geographic regions and climatic zones (Ramos Pereira et al. 2009), implying the past existence of several glacial refugia (Huang et al. 2007) including eastern and southwestern regions in China (Yan et al. 2007).
In this study, to clarify the taxonomy, distributions, and genetic differentiation of Chinese bent-winged bats, M. fuliginosus and M. magnater, we collected 125 samples across the entire range of the two species in China as well as one site in Vietnam. Firstly, we used DNA barcoding, COI, to identify the taxonomy of M. fuliginosus and M. magnater and determine their distribution ranges. Secondly, we used sequences of the mitochondrial hypervariable control region (D-loop) and NADH dehydrogenase subunit 2 (ND2) to investigate their most recent common ancestor and divergent time. Thirdly, we inferred whether climatic oscillations in Pleistocene have affected the current distribution of bent-winged bats. Lastly, we investigated the difference of morphological characters of M. fuliginosus and M. magnater.

Sampling
To encompass their entire range in China, we collected samples and morphological data of M. fuliginosus and M. magnater in China and Vietnam from 2007 to 2010. The sampling range extended across Palearctic and Oriental regions (Fig. 2). After capturing the bats in a mist net, we measured their morphological characters with digital calipers (0.01 mm) and collected wing membranes using 3-mm biopsy punchers (Worthington and Barratt 1996). The samples were preserved in absolute ethanol and stored at À20°C. The bats were released in situ. Any bats that died unexpectedly were preserved in absolute ethanol and transported to the laboratory for skull preparation and measurement. Samples were stored at the School of Environment of Northeast Normal University, Changchun, China. All field studies followed the regulations of Wildlife Conservation of the People's Republic of China (Chairman Decree [2004] No. 24) and were approved by DNA extraction, PCR amplification, and DNA sequencing Total genomic DNA of 125 individuals from 17 sampling sites was extracted using a UNIQ-10 Column animal genomic DNA isolation kit (Sangon, Shanghai, China). To analyze the genetic differentiation of bent-winged bats, we amplified the mitochondrial D-loop region for all the 125 individuals from 17 sampling sites using the two universal primers P and E as described in Wilkinson and Chapman (1991). PCR amplifications were performed in 25 lL volumes containing 2.5 lL of 10 9 PCR buffer, 2 lL dNTP mixture (10 mmol/L), 1 lL of each primer (10 lmol/L), 1 lL template DNA, and 0.5 lL Taq polymerase (5 units/lL). Samples were subjected to an initial denaturation step of 95°C for 5 min, followed by 40 cycles of denaturation at 95°C for 1 min, annealing at 55°C for 90 sec, and extension at 72°C for 2 min, followed by a final extension step of 72°C for 7 min. All haplotype sequences were deposited in GenBank under Accession Numbers KM230117-KM230241.
Based on mtDNA D-loop trees of bent-winged bats, two major monophyletic lineages were identified, corresponding to Chinese mainland and tropic region, respectively (Figs 2 and 3A). To confirm the taxonomy of those bent-winged bats, we used primers (VF1_t1 and VR1_t1) and methods described in Borisenko et al. (2008) and Francis et al. (2010) to amplify the 657 bp segment of COI bar codes of 22 individuals, which were randomly selected from each sampling sites of Chinese mainland and tropic region (Fig. 2). The resulting sequences (KM575709-KM575714, KM575717-KM575 722, KP247536-KP247545) were compared with COI sequences of M. fuliginosus and M. magnater from other researches (HM540883-HM540890).
All samples were sequenced by Sangon in Shanghai, China. The sequences were edited and aligned using Clustal 9 1.8 (Thompson et al. 1997) followed by manual adjustments.

Mitochondrial DNA analysis
Genetic diversity parameters such as haplotype diversity (h), nucleotide diversity (p), and the number of polymorphic and phylogenetically informative sites were calculated in DnaSP 5.10 (Rozas et al. 2003). Uncorrected genetic distances between and within different species were calculated using MEGA6 (Tamura et al. 2013).
The time of the most recent common ancestor (TMRCA) of the two species based on D-loop and ND2 regions was assessed using BEAST 1.6.2 (Drummond and Rambaut 2007). Because the divergence rate for the Dloop of genus Miniopterus is unknown, we used a divergence rate of 20% per million years as applied in Nyctalus bats (Petit et al. 1999). For the mutation rate for ND2, we used 1.2 and 1.8% per million years as described for Rhinolophus ferrumequinum (Flanders et al. 2009). The best substitution models, estimated using jModeltest 0.1.1, were TPM1uf + I + Γ for the D-loop and TIM2 + I + G for ND2. Because TPM and TIM models of sequence evolution are not implemented in BEAST, we used the most similar model available. The prior parameters were determined in preliminary studies. Finally, we performed runs of 30,000,000 generations, each with a burn-in of the first 10% generations, with sampling every 1,000 steps. The results were then visualized in TRACER 1.5 (Drummond and Rambaut 2007), which was also used to examine the effective sample size (ESS) of each parameter for all ESSs >1000. Neutrality tests and mismatch distribution analyses based on the D-loop sequences were used to infer population demographic events. For two mitochondrial lineages, the population demographic events of M. fuliginosus were analyzed, but M. magnater was not analyzed due to small sample size (n = 16). Fu's Fs (Fu 1997) and Fu and Li's F* and D* (Fu and Li 1993) were calculated in DnaSP 4.0. Mismatch distributions were calculated using 1,000 bootstrap replications in Arlequin. We used goodness-of-fit tests based on the sum of squared deviations (SSD) (Schneider and Laurent 1999) and raggedness index (Rogers and Harpending 1992) to test the significance of the fit of the distribution. When an expansion model could not be rejected, we estimated the time of expansion (t) from s = 2ut, where s is calculated as the time to expansion in mutational units, and u is the mutation rate per generation for the whole sequence. The values of u are equal to lgk, where l is the mutation rate per nucleotide (see above) and k is the sequence length. The generation time (g) was estimated to be 2 years (Xu et al. 2010).

Morphology analysis
To evaluate possible differentiations of the main morphological characters between M. fuliginosus and M. magnater, we analyzed forearm length (FA) for 97 samples obtained after 2007 and analyzed two major cranial parameters, maxillary third molar (M 3 ) width and greatest length of the skull (GLS) for 19 unexpected dead bodies.
We used multidimensional scaling (MDS) to arrange FA, M 3 width and GLS, between specimen pairs in a twodimensional space. From the MDS plots, the FA, M 3 width and GLS of all specimens clustered into two distinct groups with little overlap, which was consistent with the classification of three mtDNA marker. Analysis of variance (ANOVA) was therefore used to test differences between the two species. The FA of 97 individuals was measured from 13 populations (AnH, HeN, ZJ, JX, Yu-nND, ShanX, HeNWL, SC, JXSX, FJ, YueN, GD, and HaiN). The two cranial parameters were measured in 5 individuals from Hainan and 14 individuals from the Chinese mainland (AnH, YunNX, YunNW, HuN, JX, and ZJ).

Phylogeny and genetic divergence
Phylogenetic reconstructions from neighbor-joining (NJ), maximum parsimony (MP), maximum likelihood (ML), and Bayesian inference (BI) analyses produced highly concordant trees based on mitochondrial D-loop gene (Fig. 3A) Phylogenetic reconstructions from NJ, MP, ML, and BI analyses also produced highly concordant trees based on COI bar codes (Fig. 3B). The individuals randomly selected from each population of Chinese mainland (Clade I) clustered together with M. fuliginosus, and the individuals from tropical regions (Clade II) clustered together with M. magnater (Fig. 3B). The two lineages diverged by 7% based on uncorrected genetic distances, while intralineage divergences ranged from 0.1% to 1.4% for M. fuliginosus, and 0.1% to 1.8% for M. magnater. These results indicated that our samples from Chinese mainland was M. fuliginosus, and those from tropical regions was M. magnater.
The topology of ML, MP, and NJ trees based on mitochondrial ND2 was similar to those reported by Appleton et al. (2004). Similar to the D-loop trees and COI trees, all individuals clustered into two lineages, M. fuliginosus and M. magnater (Fig. 3C). M. fuliginosus from the Chinese mainland and Japan (AY169469) formed a sister species with M. magnater from tropical regions (Fig. 3C), but M. magnater from New Guinea was located in another clade (Fig. 3C). The average genetic distance between M. fuliginosus (Clade I) and M. magnater (Clade II) was 5.3%. However, M. magnater from New Guinea was relatively genetically distant to Clade I (9.7%) and Clade II (10.3%) (Fig. 3C).

Genetic diversity and estimation of divergence and expansion time
Amplification of the D-loop from 125 individuals yielded 98 haplotypes with 98 polymorphic sites and 81 phylogenetically informative sites. The total haplotype diversity (h) was 0.99 (SD = 0.002), and the overall nucleotide diversity (p) was 0.07 (SD = 0.003). Haplotype and nucleotide diversities of each population are given in Table 1. Haplotype (h = 0.99, SD = 0.003) and nucleotide (p = 0.05, SD = 0.001) diversities of M. fuliginosus were higher than those of M. magnater (h = 0.98, SD = 0.028; p = 0.04, SD = 0.003), which might be due to smaller sample sizes of M. magnater than M. fuliginosus.
TMRCA for each species was calculated using each gene. In neutrality tests of M. fuliginosus, Fu and Li's F* and D* were not significant, whereas Fu's Fs was significant; these results suggest a population expansion model for its demographic history. At the same time, Harpending's raggedness index was small (r = 0.002, P = 0.97), indicating a rapid population expansion. The inferred expansion time was 0.14 Ma (0.09-0.19 Ma).

Morphological divergence
We set up a two-dimensional MDS model to discriminate the difference of FA, GLS, and M 3 width variation between M. fuliginosus and M. magnater. The first two dimensions extracted from the MDS model described nearly all the variation in FA, GLS, and M 3 width between the two lineages (r 2 = 0.99 for FA; r 2 = 0.99 for GLS; r 2 = 0.99 for M 3 width) (Fig. 4), which revealed two distinct geographic groups, M. fuliginosus and M. magnater, similar to the phylogenetic clusters. The mean values of FA length (47.92 AE 1.06 mm; n = 81), GLS (16.09 AE 0.36 mm; n = 14), and M 3 width (6.71 AE 0.25 mm; n = 14) in M. fuliginosus were significantly slightly smaller than in M. magnater (FA: 49.50 AE 0.97 mm, n = 16; GLS: 16.91 AE 0.16 mm, n = 5; M 3 width: 7.12 AE 0.10 mm, n = 5) (ANOVA, all P < 0.01), even with little morphological overlap between the two species (Fig. 4).

Taxonomy of bent-winged bats
Our analyses indicate that DNA bar codes are an effective tool for differentiating and identifying species of bent-winged bats in China and Vietnam. Two divergent lineages (Fig. 3B) and extremely low COI genetic divergences (< 2%) between our samples and M. fuliginosus (HM540883-HM540886) or M. magnater (HM540887-HM540890) suggested that the bent-winged bats from Clade I are M. fuliginosus, and those from Clade II are M. magnater. This result is consistent with other previously studies on bent-winged bats' taxonomy (Appleton et al. 2004;Tian et al. 2004), suggesting two bent-winged bat species exist in China. However, M. magnater from New Guinea was relatively genetically distant to M. fuliginosus (9.7%) and M. magnater (10.3%) from China and Vietnam based on ND2 gene (Fig. 3C), which also suggested that further investigation is required in New Guinea to determine the correct name of bent-winged bats (Appleton et al. 2004).

Potential mechanisms of spatial distribution
In this study, these two sister bent-winged bat species correspond to different spatial distributions, the Chinese mainland for M. fuliginosus and tropical regions for M. magnater (Figs 2 and 3). Several mechanisms may be hypothesized to explain the current distribution pattern of M. fuliginosus and M. magnater, including palaeoclimatic changes (Jablonski and Whitfort 1999), climate differences (Miller-Butterworth et al. 2003;Bilgin et al. 2008), ecological attributes (Lin et al. 2014), and geographic isolation (Smissen et al. 2013).
A large body of evidence indicates that Pleistocene glacial cycles may have influenced differentiation, expansion, and genetic structure of different species (Taberlet et al. 1998;Hewitt 2000). In our study, the most recent  common ancestor of M. fuliginosus and M. magnater was dated back to 2.26 Ma on the basis of the D-loop data and to 2.37-1.69 Ma according to the ND2 data. During this time period, China had just concluded an early stage of Quaternary glaciation (2.50 Ma) and was experiencing major climatic oscillations with a dominant 0.1-millionyear cycle (Ruddiman et al. 1988;Liu et al. 2001). Climatic changes and temperature decline may have forced their ancestor to migrate southward and evolve in different refugia. Our study results indicate that the first divergence between these two species occurred either 2.26-0.58 Ma (D-loop data) or 2.37-0.42 Ma (ND2 data). This time frame extends across several glacial-interglacial stages, such as the Poyang glacial stage (1.8 Ma), the Dagu glacial stage (1.1 Ma) (Jing and Liu 1999) Pleistocene climate changes also affected population expansion events. Population expansion tests indicated that population expansion for M. fuliginosus occurred around 0.14 Ma, corresponding to the early stage (0.07-0.15 Ma) of the last interglacial period when the rising temperatures promoted population growth. Several subsequent glacial-interglacial cycles may have led to various periods of isolation and contraction of M. fuliginosus in China, consistent with its high genetic diversity (Table 1).
With respect to each species' geographic distribution, the range of M. fuliginosus corresponds to subtropical and temperate zones with a relatively cold, dry climate, whereas M. magnater has a distribution restricted to tropical coastal areas that are milder and more humid. Differing precipitation and temperature regimes in the two regions might influence the distribution of vegetation as well as insect density and composition. In our study, the two species experiencing different climatic conditions may have originally occupied different habitats and climate zones because of local ecological adaptation, and then diverged after Quaternary glaciation. This type of association between different biomes and climatic conditions has been found in two other congeneric species, M. natalensis in South Africa (Miller-Butterworth et al. 2003) and M. schreibersii pallidus in southeastern Turkey (Bilgin et al. 2008), Miniopterus manavi in Madagascar and Comoros (Goodman et al. 2009), Miniopterus fraterculus (Goodman et al. 2007), and other Miniopterus species (Christidis et al. 2014) in Madagascar, indicating that environmental conditions can influence the distribution of bats.
The ecological attributes of M. fuliginosus and M. magnater have also played an important role in modulating contemporary geographic distribution pattern. These two species have high wing loading and a high wing aspect ratio. Normally, species with this wing morphology have long-range migratory and dispersal abilities, but these two species typically roosts in caves and feeds in forests (Han et al. 2008;Hu et al. 2011). The availability of cave habitat and food resources may have caused the two species to become isolated in refugia during past ice ages. However, potential biogeographic barriers constraining other bat species (Flanders et al. 2011), such as the Qinling Mountains and the Huaihe River separating Palearctic and Oriental regions (Xu et al. 2007), did not constitute effective barriers for M. fuliginosus and M. magnater in China.

Morphology difference
Although M. fuliginosus and M. magnater can be identified based on genetic data, it is difficult to distinguish them in the field from living individuals because of their overlapping morphological characters, which is very common in Genus Miniopterus bat species, such as M. schreibersii and M. maghrebensis, or M. pallidus (Furman et al. 2009;Bilgin et al. 2012, Puechmaille et al. 2014. In this study, we found that the FA, M 3 width, and GLS were significantly different between M. fuliginosus and M. magnater, with individuals from M. magnater significantly larger than those of M. fuliginosus. This result suggested that these three parameters could be used to identify taxonomy, but still need to combine the DNA sequence because of partial overlapping. However, the width of M 3 across all individuals in our study was less than 7.3 mm (M. magnater is 7.03-7.29 mm, M. fuliginosus is 6.33-7.16 mm), whereas the width of M 3 in M. magnater was reported greater than 7.4 mm by literature (Smith and Xie 2008).

Conclusions
M. fuliginosus and M. magnater as sister species exist in China, extending across the Tropic of Cancer display a north-south distribution pattern corresponding to subtropical and temperate zones and tropical coastal areas, respectively. The TMRCA of M. fuliginosus and M. magnater could date back to the early Quaternary glacial period, with subsequent evolution occurring in different refugia. Both climate changes and their ecological attributes might have also played important roles in modulating geographic distribution pattern. Three main morphological characters, FA, M 3 width, and GLS were significantly different between M. fuliginosus and M. magnater, and the latter was significantly larger than the former.