|
nif:isString
|
-
Sample collection was carried out in collaboration with researchers located at the Thailand Institute of Scientific and Technological Research (TISTR) and directors of Wildlife Research Stations (under the administration of the Department of National Park, Wildlife and Plant Conservation in Thailand) having the required permits. This study is part of the “CERoPath project” (Community Ecology of Rodents and their Pathogens in South-East Asia: effects of biodiversity changes and implications in health ecology), ANR Biodiversity ANR07 BDIV012. CERoPath protocols (available at http://www.ceropath.org/references/rodent_protocols_book) were evaluated by the ethics committee of Mahidol University in Bangkok.
A total of 225 samples of L. neilli collected in 28 localities (Figure 1, Table S1) during a two-year survey of 122 limestone karsts distributed in all karstic regions of Thailand were used in this study (Figure S1). These 225 samples include 115 specimens from a previous study [25] and 110 new individuals, among which 50 come from eight localities (CHAI1, CHAI2, CHR, KK1, KK2, PET, SARA5, UT) unsampled in Latinne et al. [25]. Field identifications were made based on morphological criteria according to [17], [29]. Animals were live-trapped and released after taking a tissue biopsy from the ear. Skin samples were stored in 96% ethanol. Genomic DNA was extracted from skin samples using the DNeasy Tissue Kit (Qiagen Inc.) following the manufacturer's protocol.
Figure data removed from full text. Figure identifier and caption: 10.1371/journal.pone.0047670.g001 Sampling localities of the 225 L. neilli analyzed in this study. For phylogenetic analysis, six samples of L. edwardsi and L. sabanus (5 voucher specimens from the CERoPath tissue collection (R3033 (HM217400, HM217531, HQ454288, JQ081354), R3111 (HM217404, HM217534), R4222 (HM217444, HM217571, HQ454297, JQ081352), R4296 (HM217450, HM217577, HQ454292, JQ081353), R4370 (HM217451, HM217578, HQ454293, JQ081351)) [24] and one specimen collected by our team (L288 (JQ081371, JQ081319, JQ081342, JQ081355)) were chosen as outgroups (Table S2).
Mitochondrial and nuclear DNA amplification: Two mitochondrial markers, the cytochrome b gene (cytb) and the cytochrome oxidase I gene (COI), were used in this study in addition to two nuclear fragments; the intron 7 region of the β-fibrinogen gene (bfibr) and the intron 1 region of the X-linked glucose-6-phosphate dehydrogenase gene (G6pd). Sequences of cytb (n = 115), COI (n = 115) and bfibr (n = 65) were available from a previous study [25] and were associated to newly amplified sequences. The final dataset included 225 sequences of mitochondrial genes (cytb, COI) and 104 sequences of nuclear introns (bfibr, G6PD) from a subset of 104 individuals representative of the main sampling localities. Primers sets used to amplify the cytb, COI, bfibr and G6pd genes and their annealing temperatures are listed in Table S3. PCRs were carried out in 50 µl volume containing 12.5 µl of each 2 µM primers, 1 µl of 10 mM dNTP, 10 µl of 5× reaction buffer (Promega), 0.2 µl of 5 U/µl Promega Taq DNA Polymerase and approximately 30 ng of DNA extract. Amplifications were performed in thermal cycler VWR Unocycler using one activation step (94°C for 4 min) followed by 40 cycles (94°C for 30 sec, 50–58.7°C (depending on the primers, Table S3) for 60 sec, and finally 72°C for 90 sec) with a final extension step at 72°C for 10 min. Sequencing reactions were performed by Macrogen Inc. (Seoul, Korea) using sequencing on an ABI 3730 automatic sequencer. All haplotype/allele sequences have been deposited in Genbank database (GenBank Accession Numbers: cytb: HM219591–612, JQ081356–370; COI: HM219573–588, JQ081310–318; bfibr: HM219556–560, HM219562–564, HM219568, JQ081320–341; G6pd: JQ081343–350. See Table S2 for more details). For the nuclear genes, heterozygous states were identified as strong double peaks of similar height in both forward and reverse strands, or when the particular base corresponding to the dominant peak alternated on the two chromatograms [30].
The 225 samples used in this study were genotyped at twelve variable microsatellite loci using the multiplex sets and PCR conditions reported in [31]. Amplified DNA was analyzed for length variations on an ABI 3700 sequencer using GENEMAPPER 4.0 (Applied Biosystems).
Mitochondrial and nuclear DNA analysis: Sequences were aligned in BIOEDIT 7.0.9.0 [32] using the ClustalW algorithm. Haplotypes were identified using ARLEQUIN 3.11 [33]. For nuclear sequences, two alignments were created: first, heterozygous sites were coded using the IUPAC ambiguity codes (unphased dataset); and second, gametic phases of nuclear sequences were inferred using the Bayesian algorithm implemented in PHASE 2.1 [34] (phased dataset). Two runs were conducted for 1×104 iterations with the default values. Recombination in our nuclear genes was tested using the PHI test implemented in SPLITSTREE 4.11.3 [35] and the tree-based SBP and GARD methods [36] implemented online via the Datamonkey webserver [37].
Phylogenetic reconstructions were performed using the maximum likelihood (ML) and Bayesian inference (BI) approaches. Analyses were first run independently on mitochondrial and nuclear loci and then on combined datasets (cytb/COI, bfibr/G6pd (both phased and unphased datasets) and cytb/COI/bfibr/G6pd). The most suitable model of DNA substitution for each locus and dataset was determined using MODELTEST 3.0 [38] according to the Akaike Information Criterion (AIC). PhyML 3.0 [39] was used to perform ML analyses with default parameters as starting values. Robustness of the tree was assessed by 1000 bootstrap replicates. Bayesian analyses were performed with MRBAYES 3.1.1. [40]. Metropolis-coupled Markov chain Monte Carlo (MCMC) sampling was performed with 5 chains run for 3×106 generations with one tree sampled every 1000 generations. All trees obtained before the Markov chain reached stationary distribution (empirically determined by checking of likelihood values) were discarded as burn-in values. A 50% Majority-rule consensus tree was generated in PAUP 4.0b10 [41]. Median joining networks were performed with NETWORK 4.5.1.6 [42] to explore relationships among haplotypes for the mitochondrial (cytb/COI) and nuclear (bfibr/G6pd) datasets.
Genetic diversity and population differentiation: Haplotype (h) and nucleotide (π) diversities of the main lineages identified by phylogenetic analysis were estimated for the four loci independently using ARLEQUIN. Tables of nuclear allele frequency were computed with GENEPOP 4.0.11 [43] and frequency differences (genic differentiation) were tested for each nuclear locus and across all nuclear loci for all pairs of lineages containing more than three individuals with GENEPOP. The net genetic distance between lineages was computed in MEGA 4.1 [44] under the Kimura two parameters (K2P model) for the cytb dataset. Genetic differentiation among lineages containing more than three individuals was quantified by computing pairwise differentiation ΦST for the mitochondrial dataset and FST for the nuclear dataset using ARLEQUIN. A Mantel's test performed in ARLEQUIN on mitochondrial and nuclear datasets independently was used to assess the hypothesis of isolation by distance (IBD) among the four groups of populations by comparing pairwise geographic distance (log transformed) with pairwise (FST/[1−FST]).
Divergence times of L. neilli and L. edwardsi, and of the main lineages of L. neilli (approximation of the time to the most recent common ancestor-TMRCA, [45]), were estimated using Bayesian inference, as implemented in the program BEAST 1.6.1 [46] on the cytb dataset. Three fossil-based calibration points were used: (i) the split between the tribe Phloemyini and the other tribes of Murinae at 12.3 Myr (we considered Progonomys [47] as the MRCA of extant murines and assigned the oldest record of Progonomys at 12.3 Myr to the split between the tribe Phloemyini and the other tribes of Murinae rather than the Mus/Rattus split [48], [49]) (ii) the divergence between Apodemus mystacinus and all the species of the subgenus sylvaemus (A. flavicollis and A. sylvaticus) at 7 Myr [50], [51]; and (iii) the Otomyini/Arvicanthi split at 5.7 Myr [52]. Sequences of Batomys granti (AY324458) (Phloemyini) and Phloemys cumingi (DQ191484) (Phloemyini), Apodemus mystacinus (AF159394), Apodemus sylvaticus (AB033695), Apodemus flavicollis (AB032853), Otomys angoniensis (AM408343) (Otomyini) and Arvicanthis niloticus (AF004569) (Arvicanthi) were added to our dataset to calibrate the tree and constrain the age of some nodes. Analyses were performed under the TN93+G substitution model (previously estimated by MODELTEST), a strict molecular clock and a constant size coalescent population model. These priors were selected because they better fit the data than any other molecular clock and population models according to the Bayes factors calculated to compare the models [53]. However, all tested models (strict or relaxed molecular clock combined to Bayesian Skyline coalescent population model, exponential growth coalescent population model or Yule process speciation model) produced similar divergence time estimates. Three independent runs with MCMC chain length of 1.5×108 were performed, sampling every 1×104 generations. Convergence of the chains to the stationary distribution was checked using TRACER 1.5 [54]. All BEAST computations were performed on the computational resource Bioportal at the University of Oslo (http://www.bioportal.uio.no).
Genetic diversity was assessed by calculating expected (He) and observed (Ho) heterozygosities with ARLEQUIN and confirmation of Hardy-Weinberg equilibrium (HWE) was tested using GENEPOP for each locus separately and over all loci for each sampling site. The allelic richness (AR) was calculated using the rarefaction procedure implemented in FSTAT 2.9.3.2 [55]. Multi-locus Fis was calculated for each population and adjusted for multiple tests using a Bonferroni's correction with FSTAT. The proportion of null alleles (NA) at each locus and for each population was estimated with FREENA [56] and genotypes were corrected using MICRO-CHECKER 2.2.3 [57]. Tests for linkage disequilibrium between loci for each sampling site were performed with GENEPOP.
STRUCTURE 2.3.1 [58] was used to infer the number of populations (K) and assign individuals to genetic clusters independently of spatial sampling. Ten iterations were run for each value of K from 1 to 28 using an admixture model with a burn-in of 1×105 and MCMC values of 1×106. We used CLUMPP 1.1 [59] to average the results of multiple iterations for a given K. GENELAND 3.3.0 [60] was used to perform a spatial genetic analysis by integrating geographic and genetic information and to determine the most probable K. Ten runs of 1×105 MCMC iterations were performed for K from 1 to 28 with a spatial coordinates uncertainty of 1 km. The posterior probability of population membership for each individual was determined using a burn-in of 3×104 iterations. A visual output of the STRUCTURE and GENELAND results was generated using DISTRUCT [61]. Genetic differentiation among clusters containing more than three individuals was quantified by computing pairwise FST using ARLEQUIN. IBD among the four groups of populations and within groups including at least four clusters (northeast and west) was tested by comparing pairwise geographic distance (log transformed) with pairwise (FST/[1−FST]) and (RST/[1−RST]) using ARLEQUIN and SPAGeDI 1.3 [62], respectively.
We used VIP (Vicariance Inference Program) [63] to localize the main isolating barriers within L. neilli distribution range. VIP uses a phylogenetic tree (ML tree of the combined dataset) and direct geographic information data to search for the disjunctions/barriers among sister groups.
To infer ancestral areas of L. neilli populations and investigate the relative role of vicariance and dispersal in the evolutionary history of this species, we used a dispersal–extinction–cladogenesis (DEC) model of range evolution implemented in LAGRANGE [64]. The DEC model describes ancestor-descendant transitions between geographic ranges by processes of dispersal (range expansion), local extinction (range contraction), and cladogenesis (range subdivision/inheritance) and estimates likelihoods of ancestral states (range inheritance scenarios) [65]. This analysis was carried out on the chronogram inferred with BEAST limited to the ten lineages with four geographic areas (west, centre, north, northeast). We first performed an analysis without constraints to infer ancestral areas of each node of the tree. Then, to optimize the results obtained for the ancestral node of L. neilli, we constrained successively the area of this node to each combination of one, two, three or four areas (except for the combination centre+north because these two areas are not adjacent) and compared the global likelihood scores obtained for each area or combination.
Approximate Bayesian computation implemented in the software DIYABC 1.0.4.45beta [66] was used to infer the evolutionary history of L. neilli combining our mitochondrial and microsatellite datasets. Several biogeographic scenarios were compared to test whether the current groups of populations originated from the fragmentation of a common ancestral population or from one of the current groups and whether admixture events occurred among populations. According to the phylogenetic trees and Bayesian clustering analysis, four groups of populations were used (west, centre, north, northeast) and the origin of each group was tested independently in a four-step analysis using 40 scenarios (Table 1). For each group, a set of ten scenarios was designed to test its origin (see Figure 2 for a schematic representation of these ten scenarios for the NE group, similar scenarios were used for the three other groups): either from a hypothetical ancestral population (AP) (scenario NE.1), from one of the current groups of populations (scenarios NE.2 to NE.4), or a mix of two populations (admixture, scenarios NE.5 to NE.10). Range and distribution of priors for parameters used to describe these scenarios (effective population size, time of splitting or merging events, and rates of admixture in the case of merging events) are found in Table 2. For the microsatellite dataset, all one-sample summary statistics available were used in addition to FST and classification index as two-sample statistics to compare observed and simulated datasets. For the mitochondrial dataset, we used all one-sample statistics and the number of segregating sites, the mean of pairwise differences (W) and FST as two-sample summary statistics.
Figure data removed from full text. Figure identifier and caption: 10.1371/journal.pone.0047670.g002 Schematic representations of ten competing scenarios designed to test the origin of the NE group by Approximate Bayesian Computation (ABC) analysis.Ni corresponds to effective population sizes of each group. The time of events (Ti), in numbers of generation, are not sorted by increasing times but are presented randomly as their prior distribution overlapped (Table 2). Several conditions were considered: TNNE
|
![[This page as RDF]](http://data.gesis.org/somesci/static/rdf-icon.gif){kind=icon}