Contrasting Genetic Structure in Two Co-Distributed Species of Old World Fruit Bat

The fulvous fruit bat (Rousettus leschenaulti) and the greater short-nosed fruit bat (Cynopterus sphinx) are two abundant and widely co-distributed Old World fruit bats in Southeast and East Asia. The former species forms large colonies in caves while the latter roots in small groups in trees. To test whether these differences in social organization and roosting ecology are associated with contrasting patterns of gene flow, we used mtDNA and nuclear loci to characterize population genetic subdivision and phylogeographic histories in both species sampled from China, Vietnam and India. Our analyses from R. leschenaulti using both types of marker revealed little evidence of genetic structure across the study region. On the other hand, C. sphinx showed significant genetic mtDNA differentiation between the samples from India compared with China and Vietnam, as well as greater structuring of microsatellite genotypes within China. Demographic analyses indicated signatures of past rapid population expansion in both taxa, with more recent demographic growth in C. sphinx. Therefore, the relative genetic homogeneity in R. leschenaulti is unlikely to reflect past events. Instead we suggest that the absence of substructure in R. leschenaulti is a consequence of higher levels of gene flow among colonies, and that greater vagility in this species is an adaptation associated with cave roosting.


Introduction
Comparative studies taxa offer a powerful approach to elucidate how population processes and historical events have acted in shaping organisms' current distribution and genetic structure. Spatially matched taxa are often expected to show similar signatures of structure based on their shared histories [1]. For example, several studies have reported concordant patterns of genetic structure in Europe, increasing our understanding of the impact of current and past barriers to gene flow [2,3,4,5]. Conversely, disparities in genetic structure between co-distributed species can highlight differential responses to these processes often resulting from contrasting life history and ecological traits.
Dispersal ability is a key demographic force shaping natural populations [6]. Contrasting dispersal strategies, which are often linked with social structure, impact on the balance between gene flow, genetic drift, mutation and natural selection [7]. In general, species with limited dispersal abilities show more population genetic structure than do species with a tendency towards greater dispersal [8]. In bats, social structure is largely determined by roosting ecology. Natural caves provide one type of roosting habitat, but their patchy distribution in the landscape and potential to accommodate large numbers of individuals means that the number of bats using one cave can range from hundreds to thousands. Consequently, competition for local resources is expected to be high and many cave roosting bat species typically exhibit high vagility, thus allowing them to commute large distances to foraging sites [9]. Conversely, tree cavity/foliage roosting species of bats tend to live in much smaller numbers in a ubiquitous resource. This might have led to less evolutionary pressures for high dispersal. These expectations are supported by empirical data that show cavity/foliage roosting species are more susceptible to fragmentation than cave-roosting species [10].
The fulvous fruit bat (Rousettus leschenaulti) and the greater shortnosed fruit bat (Cynopterus sphinx) are co-distributed throughout much of their range in Southern Asia [11,12]. Although both species are frugivorous they exhibit markedly different roosting behavior with R. leschenaulti predominantly roosting in caves consisting of mixed sex colonies that can consist of up to 10,000 individuals [13]. In contrast, C. sphinx typically roosts in palm trees during the breeding season, and is characterized by a single male with a harem of females [7].
By occupying a similar range but having different roosting preferences, R. leschenaulti and C. sphinx provide us with an opportunity to test whether expected differences in dispersal show a correlation to population genetic structure. In this study, we aim to compare the population genetic structure of R. leschenaulti and C. sphinx over a wide geographical area using both mitochondrial (cytochrome b gene) and microsatellite markers. We predict that R. leschenaulti will have a high dispersal ability that will be characterized by the identification of low population genetic structure across its range, while the opposite will be found for C. sphinx due to its predicted limited dispersal ability.

Ethics statement
Sampling was approved by the Administrative Panel on Laboratory Animal Care (approval number 2009001) of the Guangdong Entomological Institute in China, which also incorporates the South China Institute of Endangered Animals.

Sample collection
Rousettus leschenaulti was sampled at four localities in China (Maoming, Wuming, Haikou, Menglun) and one in India (Cheranmahadevi). Cynopterus sphinx was sampled at six localities in China (Guangzhou, Zhongshan, Jiangmen, Beihai, Haikou, Xishuangbanna), one in India (Mandiyoor) and one in Vietnam (Pu Huong) (Table 1, Figure 1). To avoid over-representation of potential relatives, bats of both species were captured using mist nets at foraging grounds (or for Rousettus, also near to roost entrances) rather than by hand-netting. Tissue samples were taken from the wing-membrane using a 3-mm diameter biopsy punch and stored in 80% ethanol until processing [5]. All bats were released at the site of capture DNA extraction and microsatellite genotyping DNA was extracted using DNeasy Tissue kits (Qiagen). R. leschenaulti was genotyped at eight microsatellite loci while C. sphinx was genotyped at six (see Table S1). To amplify microsatellite loci, polymerase chain reactions (PCRs) were undertaken using total reaction volumes of 10 ml, each containing 50-100 ng genomic DNA, 0.25 mM of each primer, and 1x PCR buffer containing 2 mM MgCl 2 , 0.2 mM of each dNTP, and 0.25 U hot-start Taq DNA polymerase (Qiagen). Forward primers were 59-fluorolabelled (MedProbe). PCRs were performed using a PTC-200 thermal cycler (MJ Research) with the following profile: denaturation at 95uC for 15 min, 35 amplification cycles (94uC for 30 s, annealing temperature for 30 s, 72uC for 30 s) and an extension step at 72uC for 20 min. PCR products were genotyped using an ABI 3100 DNA sequencer (Applied Biosystems) and alleles were sized using the programs GENESCAN version 2.1 and GENO-TYPER version 2.5 (Applied Biosystems).

mtDNA sequencing
Five to eight individuals from each population (total n = 88) were amplified at the cytochrome b gene (,1100 bp) using the published primers cy1 (59-AAA TCA CCG TTG TAC TTC AAC-39) and cy2 (59-TAG AAT ATC AGC TTT GGG TG-39) [14]. PCRs were performed using a PTC-220 thermal cycler (MJ Research) in 50 ml reaction volumes, each containing 25 ml of 26 ExTaq polymerase (Takara), 0.5 ml of DNA templates (50 mg/ml), and 2 ml of each primer (10 mM). For PCR conditions see the microsatellite genotyping section above. PCR products were sequenced using Big Dye Terminator kits (Applied Biosystems) on an ABI 3730 automated sequencer with both primers.

Microsatellite statistical analyses
We tested for deviations from linkage equilibrium between loci in FSTAT version 2.9.3.2 [15], and from Hardy-Weinberg equilibrium (HWE) for each population and locus using the Markov chain method (10,000 dememorization steps, 10,000 batches and 10,000 iterations per batch) in GENEPOP version 3.3 [16]. Population values of the mean number of alleles per locus and heterozygosity (observed and expected) was calculated using GENEPOP. We also estimated allelic richness (R S ) per locus per population, which accounts for unequal sample sizes, in FSTAT. To estimate the power of the markers, we calculated the total probability of identity' (PID) for each locus and for all loci combined in GENALEX [17].
To quantify genetic structure, we calculated pairwise F ST values using the software GENETIX version 4.02 [18]. To test for a phylogeographic signal in the data, we also derived pairwise R ST values, which accounts for the stepwise mutation model of microsatellite alleles [19], and compared these to R ST values in which allele sizes were permuted (1000 times) among alleles ( p R ST ) [20] using the software SPAGeDi version 1.3 [21]. The extent to which R ST exceeds p R ST is thus a measure of the extent to which phylogenetic distance among ordered alleles explains the observed pattern of genetic differentiation. We then assessed the relative contribution of phylogeography between the two species by using a paired sample t-test on equivalent pairwise values of R STp R ST .
We also plotted linearized pairwise values of both F ST and R ST against log-transformed geographic distances (km) to test for isolation by distance (IBD), using a Mantel test with permutations (1000 times) to assess significance in GENEPOP. An analysis of molecular variance (AMOVA) was used to examine genetic variation within and among the different populations, based on all samples, with and without the Indian colony (AMOVA I and II, respectively). We assessed whether the derived indices were significantly different from zero using a permutation procedure (5,000 iterations) in the software ARLEQUIN version 3.1 [22].
Population structure was further investigated using Bayesian clustering analysis. We estimated the likelihood of different numbers of clusters (K) in the data using the method in STRUCTURE version 2.0 [23]. Five independent runs (burn-in of one million and one million MCMC steps), were undertaken for each value of K from 1 to 8. We applied the admixture symmetric similarity coefficients (SSC) among replicate runs within each value of K [24], using the Greedy algorithm of CLUMPP version 1.1.1 [25] in which groups of runs with a SSC$0.8 were identified and combined. Summary outputs for each value of K were then displayed graphically using the software DISTRUCT version 1.1 [26].

Sequence analysis
Cytochrome b sequences were aligned in BioEdit version 7.0.5 [27] and both haplotype diversity (h) and nucleotide diversity (p) were calculated for each population in the software DNASP version 4.0 [28]. For each species, the relationship among haplotypes was assessed using a network method based on statistical parsimony. Haplotype networks were constructed in the program TCS version 1.21 [29] with the maximum number of mutational differences justified by the parsimony limit of 0.95.
Tests for genetic differentiation among samples (n$5), calculated using pairwise W ST values and tested for significance by permutation (10 000). AMOVAs were used to examine genetic variation following the methods described above. Finally, to test whether any detected differences in population structure between the two species could be explained by contrasting demographic histories, we performed sequence mismatch distribution analyses in ARLEQUIN. Mismatch distributions are typically ragged or multimodal for populations at stationary demographic equilibrium, but smooth or unimodal for populations that have undergone a demographic expansion [30]. Goodness of fit tests for a model of population expansion were calculated from the sum of squared deviation (SSD) and the raggedness index (r), and significance was assessed by bootstrapping (10 000 replicates). Where evidence of population expansion was found, the expansion time in generations (t) was derived following t = T/2u, where T (tau) is a parameter of the time to expansion in units of mutations, and u is the mutation rate per generation for the DNA sequence. We used a mutation rate of 2% per Myr [31] with a generation time of 2 years, based on age of first breeding for most insectivorous bat species [32].

Genetic diversity
Samples of both Rousettus leschenaulti and Cynopterus sphinx showed no evidence of linkage disequilibrium between loci, and there was no significant deviation from Hardy-Weinberg equilibrium detected in either species following Bonferroni correction for multiple tests.

Genetic population structure
Levels of population genetic differentiation appeared to differ between the two species. Global F ST for C. sphinx was considerablly greater (F ST = 0.0235, P,0.001) than that recorded for R. leschenaulti (F ST = 0.0067, non-significant). Pairwise F ST values among C. sphinx samples ranged from 0.0090 (non-significant) for Pu Huong vs. Jiangmen to 0.0801 (P,0.001) for Jiangmen vs. Mandiyoor with 22 out of 28 pairwise values found to be significant (Table S2). Pairwise values among R. leschenaulti samples were significantly different in just three out of ten comparisons (all involving comparisons with Cheranmahadevi and Haikou) (Table  S3). Although C. sphinx was sampled at more localities than was R. leschenaulti, thus precluding a meaningful comparison of their overall level of differentiation, the above species trend was also evident for those comparisons that were exactly or approximately geographically matched.
Global R ST was not significantly different from global p R ST (P,0.05) in either species indicating that stepwise mutation had not contributed to observed differentiation and that F ST is a more suitable estimator in this study. Hierarchical AMOVAs based on F ST revealed that more genetic variance was partitioned among individuals within populations than among populations in C. sphinx (96.62% versus 3.38% respectively) and R. leschenaulti (99.48% versus 0.52% respectively). Removing the populations from India did not significantly alter the result. C. sphinx was also found to exhibit a significant pattern of isolation by distance based on pairwise F ST (R 2 = 0.778; P,0.01). This analysis was not undertaken for R. leschenaulti due to the smaller number of pairwise comparisons.
In agreement with F-statistics, Bayesian clustering analyses also revealed differences between the two taxa ( Figure 2). For C. sphinx, K = 4 showed highest probability, however, forcing values of K from 2 to 8 recovered hierarchical population structure. At K = 2, Mandiyoor in India and Xishuangbana in SW China formed a single cluster with all other populations forming a second cluster, though partial admixture was observed. At K = 3, Mandiyoor and Xishuangbana remained as a single cluster, while the five populations from China showed partial membership to the remaining two clusters, one of which being dominated by the colony from Haikou (Hainan Island). Pu Huong appears to share mixed ancestry from all the sampled populations. At K = 4 the Mandiyoor and Xishuangbana populations split with Xishuangbana largely forming a new cluster, although some admixture is seen, while the populations from China and Pu Huong show very little structure and show partial membership to all the clusters apart from the one for Mandiyoor. Clustering analyses revealed no

MtDNA sequence analysis
We sequenced 1100 bp of cytochrome b in 29 individuals of R. leschenaulti and recorded 27 haplotypes and 52 polymorphic sites (GenBank accessions FJ489931 -FJ489957). For C. sphinx we sequenced 1075 bp of cytochrome b in 59 individuals and found 34 haplotypes with 71 polymorphic sites (GenBank accessions FJ489958 -FJ489992). No insertions or deletions were observed in either species. For each sample, the number of haplotypes, haplotype diversity, the number of polymorphic sites, average number of differences and nucleotide diversity are given in Table 1.
Parsimony-based haplotype networks revealed differences between the two taxa ( Figure 3). For R. leschenaulti, haplotypes formed a single network at a 95% parsimony threshold, and were scattered with respect to locality. No evidence of any structure was present.
The most common haplotype was geographically widespread, recorded in both Maoming and Cheranmahadevi. In C. sphinx, a parsimony threshold of 95% yielded two subnetworks, one comprising only haplotypes from Mandiyoor and the other comprising the remaining haplotypes. The most common haplotype was found mainly in Guangzhou, Zhongshan and Jiangmen but was also recorded in Pu Huong. Based on its interior position and high level of connectedness, it is likely to be ancestral. The lack of monophyly seen in haplotypes from Haikou, Beihai and Xishuangbanna points to past mixing among these localities.
No genetic differentiation was detected based on mtDNA data (global exact test, P.0.05) for R. leschenaulti populations. In C. sphinx, significant genetic differentiation was detected (global exact test, P,0.05) with pairwise values indicating that this is caused by the population sampled at Mandiyoor which was consistently differentiated from all other colonies sampled in China and Vietnam (P,0.05). Apart from the tests including Mandiyoor, no significant differentiation was detected between any other pairwise values. AMOVAs carried out for R. leschenaulti showed no significant genetic variance among populations, with or without Cheranmahadevi, with 3.90% of the variance explained by amongpopulation variation, and 96.10% by within-population variation (Table S4). AMOVAs for C. sphinx revealed that significant genetic variance was attributable to the two hierarchical levels examined (among and within populations). When all eight populations were included in the analysis (AMOVA I), 55.11% of the total genetic variance was explained by among population differences, and 44.89% by within population differences, however, when Mandiyoor was excluded (AMOVA II), these values were reduced to 3.04% and 96.96%, respectively (Table S4).
Mismatch distribution analysis for R. leschenaulti failed to reject the model of population expansion (P SSD .0.05 and raggedness index P R .0.05) with an estimated timing of expansion occurring around 200 000 years BP (95% CI: 130 000-250 000 years BP). In C. sphinx, separate mismatch distribution analyses undertaken for all localities also failed to reject the expansion model (P SSD .0.05 and raggedness index P R .0.05), with an estimated timing of expansion of around 107 000 years BP (95% CI: 16 000-560 000 years BP).

Discussion
In this study we applied mtDNA and microsatellite analyses to characterize the phylogeographic histories of co-distributed populations of Rousettus leschenaulti and Cynopterus sphinx. We show that in spite of their similar ranges, these two fruit bat species exhibit highly contrasting patterns of genetic structure, and that these differences are more likely to reflect differences in vagility rather than demographic history.
Populations of R. leschenaulti were characterized by a lack of genetic structure at both classes of molecular marker. Our MSN network showed that haplotypes were not geographically structured while tests of AMOVA indicated that 96% of the genetic variance was due to within-population variation. Microsatellite data also revealed complete mixing of localities based on clustering analyses of genotypes which was further supported by limited evidence of genetic differentiation from F ST statistics.
In contrast, haplotype data for C. sphinx showed a lack of genetic structure among all populations sampled across China and Vietnam but significant genetic differentiation between each of these populations and that from Mandiyoor (India). These results were also reflected in the distribution of haplotypes in the MSN network, in which the Indian haplotypes formed a monophyletic sub-network. Greater substructure was detected from microsatellite data, with significant genetic differentiation seen among most pairs of populations, including 11 pairwise comparisons within China. Bayesian clustering of genotypes also revealed clear subdivision. Specifically, we found that samples from Yunnan in Southwest China showed greater similarities with samples from India than those from elsewhere in China and Vietnam. Sampling of additional populations in the adjoining regions would help to verify these results and delimit the geographical extent of these haplotype affiliations.
A lack of genetic structure can reflect either gene flow or a recent expansion. We tested whether broad differences between our two focal taxa could be explained by differences in their demographic history. Although we detected expansions in both taxa, we estimated that the timing of population growth in the more genetically homogenous R. leschenaulti occurred around 200,000 years BP, which was earlier than the date of inferred growth for C. sphinx (107,000 years BP). Our results therefore indicate that the lack of genetic structure in R. leschenaulti is unlikely to reflect a recent expansion but instead results from a high level of gene flow among localities, both contemporary and in the past. In C. sphinx, however, the comparatively higher level of structure appears to have resulted from restricted gene flow since a more recent expansion.
Although our results are based on few and unevenly distributed populations, and so should be treated cautiously, the observed  [33,34]. Therefore, depending on whether mating occurs before or during migration this could provide further opportunities to promote gene flow [but see: [35]]. In contrast, C. sphinx forms smaller colonies, often in palm trees, and typically comprising a single male with  multiple females [36]. A previous study of this taxon at a landscape scale revealed a high variance in male reproductive success yet only moderate genetic differentiation among groups [37]. In our study, mist netting of C. sphinx at foraging sites rather than at roosts was undertaken to reduce the potential over-representation of relatives in our samples; however, we cannot rule the possibility that this has contributed to observed subdivision. In addition to the overall trend of restricted gene flow, our results of C. sphinx from Yunnan also revealed a marked discontinuity in allele frequencies within China and Vietnam, pointing to a phylogeographic signature. In particular, we found that bats sampled from Yunnan clustered with more geographically distant individuals from India rather than with other bats in China and Vietnam. In recent years, deep genetic subdivision in Southwest China, with associated suggestions of multiple glacial refugia and contact zones in this region, have been reported for several taxa, including the greater horseshoe bat [38] and the ringnecked pheasant [39]. Although the sampling regime in this study cannot resolve or explain such subdivisions, our results do hint at the possibility of multiple refugia. Our network analyses of the cytochrome b sequence data based on a parsimony threshold of 95% confirmed the monophyly of the Indian samples previously described [40]. This distinction was also reflected in the contrasting degree of population subdivision recorded when India was either included (Table S4, 55.11%, AMOVA I) or excluded (among population, 3.04%, AMOVA II) from our analyses, so supporting earlier suggestions that C. sphinx in India and China probably belong to different subspecies [41]. Interestingly, however, no such pattern was seen in the corresponding AMOVA results based on the microsatellites (among population variance 3.38% and 2.52% with and without the Indian sample, respectively). Given that both of the corresponding pairwise Fst values were not especially high, one possible explanation for this result is male-biased gene flow, either current or in the past. Alternatively, and arguably more reasonable given the scale involved, the low Fst between India and China could reflect microsatellite allele size homoplasy, which can sometimes result in erroneously low Fst values between geographically distant samples [e.g. [5]].
Comparative studies of co-distributed taxa provide a powerful means of disentangling the relative impact of biotic and abiotic factors on population genetic structure [e.g. [42,43]]. A recent study of five passerine bird species from the Tibetan Plateau revealed that species-specific demographic histories depended on habitat requirements and, therefore, local distributions [44]. In contrast, the bats in our study have broadly similar habitat requirements, and thus the observed differences in genetic structure are probably more likely to be shaped by roosting and behavioural differences that hold across a wider sampling area. Our findings support a number of papers that have linked social structure to genetic differentiation [e.g. [45]], and also corroborate suggestions that colonial cave roosting bats are likely to show more dispersal and gene flow than do tree roosting species [10].