Email updates

Keep up to date with the latest news and content from BMC Evolutionary Biology and BioMed Central.

Open Access Highly Accessed Research article

Influence of ancient glacial periods on the Andean fauna: the case of the pampas cat (Leopardus colocolo)

Daniel Cossíos1*, Mauro Lucherini2, Manuel Ruiz-García3 and Bernard Angers1

Author Affiliations

1 Département de Sciences Biologiques, Université de Montréal, C.P.: 6128, Succ. Centre-Ville, Montréal, H3C 3J7, Canada

2 GECM, Cátedra de Fisiología Animal, Departamento de Biología, Bioquímica y Farmacia, Universidad Nacional del Sur-CONICET, San Juan 670, 8000 Bahía Blanca, Argentina

3 Departamento de Biología, Facultad de Ciencias, Pontificia Universidad Javeriana, Cra 7A No 43-82, Bogotá DC, Colombia

For all author emails, please log on.

BMC Evolutionary Biology 2009, 9:68  doi:10.1186/1471-2148-9-68


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2148/9/68


Received:12 August 2008
Accepted:30 March 2009
Published:30 March 2009

© 2009 Cossíos et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

While numerous studies revealed the major role of environmental changes of the Quaternary on the evolution of biodiversity, research on the influence of that period on current South-American fauna is scarce and have usually focused on lowland regions. In this study, the genetic structure of the pampas cat (Leopardus colocolo), a widely distributed felid, was determined and linked to ancient climate fluctuations on the Andean region.

Results

Using both mitochondrial sequences and nuclear microsatellites, we inferred the existence of at least four groups of populations in the central Andes, while other three localities, with little sample sizes (n = 3), presented differences in only one of these markers. The distribution of these groups is correlated to latitude, with a central area characterized by admixture of numerous mitochondrial clades. This suggests colonization from at least three glacial refuges and a contact zone between 20 degrees and 23 degrees S following a glaciation event. The similar coalescence times of the mitochondrial haplotypes indicated that the major clades split approximately one million years ago, likely during the Pre-Pastonian glacial period (0.80 – 1.30 MYA), followed by a demographic expansion in every clade during the Aftonian interglacial period (0.45 – 0.62 MYA). Interestingly, this structure roughly corresponds to the current recognised distribution of morphological subspecies.

Conclusion

The four groups of populations identified here must be considered different management units, and we propose the three localities showing differences in only mtDNA or ncDNA as provisional management units. The results revealed the influence of ancient climate fluctuations on the evolutionary history of this species. It is expected that the other species of land vertebrates with a smaller or similar mobility have been affected in the same manner by the glacial and interglacial periods in the central Andes

Background

Previous bio- and phylogeographic studies in the northern hemisphere revealed the major role of the Quaternary in the evolution of biodiversity. Important climatic fluctuations resulted in alternations of glacial and interglacial periods that profoundly influenced the distribution of habitats, dispersal and isolation among populations [e.g. [1-4]]. In South America, climate fluctuation periods were equivalent to those known for the northern hemisphere [5,6]. During the glacial periods, glacier boundaries descended hundreds of meters, and advanced dramatically in the southern part of the continent [5,7,8]. In the central Andes, although late Quaternary glaciations eliminated the material evidence of glacier advances older than 200 000 years [5], research has shown that the distribution of plant species and habitats cycled with climatic changes, descending to lower elevations during periods of reduced temperature [9,10]. However, studies on the influence of the Quaternary on current South-American biodiversity are scarce and have usually focused on lowland regions [11,12].

Species with broad ranges of distribution can be informative regarding glacial refugia, dispersal pathways and contact zones [13,14] and, then, provide useful insights about the role of climatic changes on biodiversity. The pampas cat (Leopardus colocolo) displays a large distribution in South America, from Central Ecuador to Patagonia, in a variety of habitats [15,16]. Based on morphological characteristics, between 8 and 11 subspecies are currently recognized [15,17,18] and the scission of this taxon into three different species was proposed by García-Perea [17]: L. colocolo for populations distributed on the western slope of the southern Andes; L. pajeros distributed along the Andes; and L. braccatus found to the east of the Andes, in Brazil, Uruguay and Paraguay. Previous studies performed on mitochondrial genome [19,20] also revealed that the Andean pampas cat populations are genetically structured and may have experienced significant and lengthy periods of isolation and reduced gene flow.

In spite of its wide distribution, the pampas cat is one of the less known felids [21] and its status is affected by a variety of threats, comprising habitat loss and fragmentation, hunting for traditional reasons and decline of prey populations [22,23]. Lack of evaluation of its conservation status through its range causes the pampas cat to be considered as vulnerable [24] and to be included in the IUCN Near Threatened (NT) category [25].

The aim of the present study was to link the genetic diversity of the pampas cat throughout the central Andes to ancient climate fluctuations. To address this objective, the genetic structure of the pampas cat in 19 localities over a distance of more than 4000 km was inferred with both mitochondrial and nuclear genomes, using non-invasive sampling.

Results

Species identification and mtDNA variability

A total of 39 skin samples and 532 faecal samples from 19 localities were analyzed (Figure 1; Table 1). Of the faecal samples, 406 were unambiguously assigned to pampas cats according to the mtDNA 16S region. The remaining samples were assigned to Andean cats Leopardus jacobita (40), domestic cats (11) or canid species (32). A small number of faecal samples (43) failed to be amplified and could not be identified. All skins were assigned to pampas cats.

Table 1. Sampled localities, geographical coordinates and sample sizes

thumbnailFigure 1. Pampas cat sampled localities. The locality numbers correspond to those in Table 1. The grey and yellow areas refer to the approximate distribution of the pampas cat subspecies mentioned in this research, as proposed by García-Perea (1994).

The length of two mitochondrial control region segments combined varied between 336 and 351 bp. Two variable mononucleotide repeats (T)3–7 (C)3–11 were identified but not considered for the analyses. Without considering the RS2 region composed of a variable number of large tandem repeats, the pampas cat HVS-I region contains between 178 and 190 bp (Figure 2). The range known for other felid species is 231–245 bp [26-28].

thumbnailFigure 2. Pampas cat hypervariable sequence I (HVS-I). The position of the primers used in this study is shown. In this example, the repeated sequence 2 (RS2) contains three tandem repeats named here a, b and c. Although the primer H1rev anneals on the a, b and c repetitions and the primer H2for anneals on the a and b repetitions, only their positions on the sequences used in this research are shown.

SSCP analyses revealed a single allele per sample, excluding the presence of mtDNA sequences transferred into the nuclear genome (Numts) [28]. The SSCP survey of the 406 samples and the sequencing revealed a total of 41 HVS-I haplotypes and 94 variable sites. Phylogenetic relationships inferred among these haplotypes were consistent for both NJ and ML methods. The haplotypes clustered into four major clades, named hereafter A to D, strongly supported by bootstrap values (Figure 3). However, relationships among clades were not completely resolved.

thumbnailFigure 3. Neighbour-joining tree of the observed HVS-Isequences. Labels correspond to haplotype identification numbers. Only bootstrap values over 50% of 1000 bootstraps are shown, corresponding to the NJ/ML analyses. Two ocelots (Lpa1 and Lpa2) and two margays (Lwi1 and Lwi2) were used as outgroups. Scale bar represents an interval of Tamura-Nei genetic distance.

Sequencing of the NADH-5 and ATP-8 genes revealed 29 sites of additional variation for the 19 individuals selected from our sample. These individuals were selected randomly, with the constraints of choosing individuals with different HVS-I haplotypes and covering all the major HVS-I clades. The phylogenetic tree inferred using both the NADH-5 and ATP-8 genes was fully congruent with the one performed with HVS-I and allowed the four major clades to be recovered, while the resolution was lower. When individuals from previous studies [20,29] were included into the present dataset, a total of 41 variable sites was detected for NADH-5 and ATP-8 genes. Most of these additional individuals clustered within clade B (15) and D (3) (Figure 4) while individuals geographically distant from the Andean region, located in central Chile and Brazil, formed two additional clades. Interestingly, individuals presenting HVS-I haplotypes of the clades A and C formed two clusters not identified by previous studies.

thumbnailFigure 4. Neighbour-joining tree of the combined observed NADH-5 – ATP-8 sequences. Sequences obtained in this study are marked with a star and their labels indicate their haplotype and, between parentheses, their localities. Additional sequences from Chile, Bolivia and Brasil were previously released by Johnson et al. (1998) and Napolitano et al. (2008), and their labels correspond to individual identification numbers. For each cluster, the correspondent clades of the HVS-I are indicated for the individuals sequenced in this research, as well as the geographical origin of the samples analysed by the other authors. Only bootstrap values over 50%of 1000 bootstraps are shown, corresponding to the NJ/ML analyses. Two ocelots (Lpa11 and Lpa31) and two margays (Lwi22 and Lwi36), from Johnson et al. (1998), were used as outgroups.

Individual identification and population diversity

Microsatellite amplifications provided results for 290 (71%) of the 406 pampas cat faecal samples and for the 39 skin samples (100%) for a total of 329 samples. The microsatellites were highly variable, with 10–21 alleles per locus (Table 2) and the probability of sampling two different individuals with the same genotype ranged between 4.80 × 10-3 and 2.88 × 10-16. Unique multilocus genotypes were recorded for 99 faecal samples (30%). However, in several cases, two or more samples displayed the very same multilocus genotype indicating that the same individual had been sampled several times. According to the low probability of obtaining the same genotype in different individuals, samples with the same genotype were assigned to a unique individual, providing a final sample size of 199 pampas cat individuals. The number of individuals per locality varied from 5 to 30, except for 6 localities that had a sample size of 4 individuals or less (Table 1). Localities with sampling size lower than four individuals were not included in the following analyses, unless mentioned in the text.

Table 2. Haplotype and microsatellite diversity in pampas cat

The total number of alleles per population ranged from 17 to 48 (Table 2), but was correlated to sampling size (Pearson's r = 0.915, P < 0.0001). The total number of alleles estimated for 5 individuals per population (allelic richness, FSTAT 2.9.3) [30] ranged from 17.16 to 24.85, showing little variation between populations. The expected heterozygosity ranged from 0.35 to 0.93 (Table 2) and was not correlated to sampling size (r = 0.013, P = 0.960). None of the localities displayed deviation from HW expectations, indicating that no more than one population was sampled per locality.

For the mtDNA control region, the number of haplotypes per population ranged from 2 to 14 (Table 2) and was also correlated to sampling size (r = 0.884, P < 0.0001). Haplotype diversity values ranged from 0.60 to 0.93 between the sampled localities, and nucleotide diversity varied by one order of magnitude between 0.0059 and 0.0519.

Population structure

The best value of the ln Pr(X|K) obtained with the program STRUCTURE on microsatellite data corresponded to K = 3 population groups. These groups included the localities (2–6, 8, 9), (11, 14) and (15–18). Relationships among populations inferred from the microsatellite data supported this grouping, although bootstrap values between populations 2–6, 8 and 9 were low (Figure 5A). AMOVA analysis also supported this structure as the one displaying the highest variation among groups on either mitochondrial or microsatellite data (Table 3).

Table 3. AMOVA results for groupings of populations

thumbnailFigure 5. Population structure inferred from microsatellite data. A. Neighbour-joining phylogenetic tree of the sampled localities with five or more individuals, constructed using Dce distances and 200 bootstraps on locus. B. Membership coefficients inferred with the program STRUCTURE (Pritchard et al 2000) for K = 3. Each individual is represented by a column and each of the three inferred population groups is represented by a colour. A star indicates localities with less than five individuals, that were not included in the neighbour joining tree showed in (A).

The distribution of the groups of populations appeared to be clearly correlated to latitude (Figure 6). The first group occupies the area north to 18°S and is formed almost exclusively by individuals of the clade A. The second group, distributed between 20° and 23°S, presents a great proportion of individuals of clade C, although clades B and D are present too. The third group, distributed south to 25°S, contains principally clade D individuals, with some clade B and C individuals in its northern region. The northernmost and the southernmost sampled localities are represented by private haplotypes of the clades A and D, respectively.

thumbnailFigure 6. HVS-I haplotypes distribution and frequencies. Diameter of circles are proportional to the number of individuals. Numbers assigned to the haplotypes correspond to those on Figure 4.

Including all localities did not modify the results provided by STRUCTURE: (1–10), (11, 13–14) and (12, 15–19) (Figure 5B) and the population structure remains strongly correlated to geography except for population 12. While geography and mtDNA data suggested including this population in the group (11–13–14), the microsatellite data suggested a closer affiliation with the southernmost group (15–19).

Evolutionary and demographic history

Based on the HVS-I sequences, the coalescence times of the clades [A, B, C and D], [A – B], and [C – D] were estimated to 1.52 MY ago (with a confidence interval from 1.05 to 2.21 MYA), 1.37 MY ago (0.95 – 2) and 1.17 (0.81 – 1.7) MY ago respectively. Interestingly, the radiation of the clades A, B, C and D occurred almost simultaneously, with coalescence times of 0.43, 0.39, 0.54 and 0.44 MY ago, respectively. The results for Tamura-Nei distances between clades, coalescence times and confidence intervals are presented in Figure 7.

thumbnailFigure 7. Divergence times between pampas cat clades and their correspondence with glacial and interglacial periods. Values were estimated considering a divergence time of 2.91 MY between the pampas cat and the ocelot and the margay, with a confidence interval between 2.02 and 4.25 MY, as estimated by Johnson et al. (2006). Interglacial periods are highlighted in grey.

Discussion

Genetic diversity

Previous research has shown that pampas cats present a large genetic diversity at the species level [19]. Results of the present study indicate that a large diversity is also present at populations level and that it is comparable to other non-endangered felid species, like the ocelot and margay [27]. Although we reported a higher HVS-I genetic diversity for some pampas cat populations (P = 2.35 – 19.70; Table 2) than for those species (P = 3.71 – 14.73) [27] this can be explained by bigger samples and the presence of highly divergent clades in the most diverse pampas cat populations (Figure 6).

Genetic structure and geographic distribution

Along the Andean region, analysis of mitochondrial and nuclear genomes revealed the presence of three groups of populations. These groups show a strong geographical structure, with latitudinal separations at 18°–20°S and 23°–25°S. Including results of previous studies based strictly on mitochondrial DNA indicated the existence of a fourth group in northern Chile [20], as well as the two groups out of the Andean region, in Brazil and central Chile [19]. This genetic structure roughly corresponds to the distribution of morphological subspecies described by García Perea [17]. Two of the museum samples, both from locality number 8, were identified as L. c. garleppi using the descriptions done by García-Perea [17], although it was not possible to directly assign a genotype to a morphotype for other localities due to the nature of the samples. However, according to the type localities of the described subspecies, the pampas cat groups can be equivalent to garleppi (Clade A, 9°–18°S), budini (20°–23°S), pajeros (25°–38°S) and wolffshoni (Clade B northern Chile) subspecies (Figure 8).

thumbnailFigure 8. Distribution of the seven Management Units (MUs) proposed for the pampas cat in the Andean and Argentinean pampas regions, as defined by mtDNA and microsatellite analysis. The proportion of the mtDNA clades and the sample size are shown for each MU. MUs that are encircled with the same colour presented differences in only mitochondrial or nuclear markers, and those with only 3 individuals sampled are proposed as provisional. The MU without indicated sample size was proposed on the basis of Napolitano et al (2008) mtDNA results. The grey and yellow areas refer to the approximate distribution of the pampas cat subspecies mentioned in this research, and their names are indicated. Dots correspond to the sampled localities.

The geographical location of populations 12 and 19 coincides with the subspecies steinbachi and crucinus, respectively, while population 1 is located between the supposed distribution ranges of thomasi and garleppi (Figure 8). These three localities are characterized by private HVS-I haplotypes or ncDNA that differs from that of adjacent populations, and had samples of only 3 individuals. Since a small sample size per locality can affect the results of the STUCTURE program [31] and other analyses, the correspondence between these localities and subspecies needs to be validated by further investigations. The further assignation of these localities to subspecies can be of conservation value since many authors [i.e. [18,19]] do not recognise steinbachi and crucinus as valid taxa. As this research focused in the central Andean region, further research must be made to describe the genetic structure of the entire species.

The Andean area between 18° and 23°S, approximately, corresponds to an extremely arid belt that separates a northern area, with summer rains, from a southern one, with winter rains [32,33]. This zone is considered to be a barrier between subspecies of other land mammals such as the vicuna Vicugna vicugna [34] and the lesser grison Galictis cuja [35], and represents a distribution limit for other species, like the long-tailed weasel Mustela frenata [36]. In contrast, the genetic structure of the puma Puma concolor from the high Andes is not affected by this barrier [37], probably due to a higher dispersal capacity. The pampas cat populations distributed between 18°S and 25°S display three different phenotypes [17] and an important admixture of the clades typical of adjacent areas, suggesting the arid belt as a contact zone.

Influence of Pleistocene

The topology of the mtDNA trees enabled us to infer two periods of prime importance for the demographic history of the pampas cat. The lack of resolution for the relationships between pampas cat clades (as well as among haplotypes within clades) suggests a rapid radiation.

A first episode comprised the split of the clades A, B, C and D and occurred between the end of the Bramertonian Interglacial (1.30 – 1.55 MYA) and the beginning of the Pre-Pastonian glacial period (0.80 – 1.30 MYA; Figure 7). The Pre-Pastonian corresponds with the most extensive glaciations in southern South America [38]. This period also coincides with the estimated date of the divergence between the major clades of the Andean bird genus Muscisaxicola [39], suggesting this period as an important phase in the diversification of the Andean fauna.

The very similar coalescence times of the current haplotypes of the clades A, B, C and D suggest these events took place simultaneously, during a second demographic episode. These splits were likely to be the result of demographic expansions caused by geographically extended phenomena associated with climate change. The calculated mean time for these events overlapped the end of the Kansan glacial period (0.30 – 0.45 MYA) and the Aftonian interglacial (0.45 – 0.62 MYA), the interglacial period being the more likely moment of occurrence for these demographic events (Figure 7).

Interestingly, the divergence time between clades A, B, C and D is similar or longer than between some felid species, such as the Iberian lynx Lynx pardinus, Canadian lynx Lynx canadensis and Eurasian lynx Lynx lynx (between 1.18 and 1.61 MY), kodkod Leopardus guigna, little spotted cat Leopardus tigrinus andGeoffroy's cat Leopardus geoffroyi (0.74–0.93 MY) and ocelot and margay (1.58 MY) [40].

High diversity within clades, times of divergence and monophyletic origin in most of the regions indicated long-term isolation into distinct refuges (glacial)/regions (interglacial) rather than dispersal from a unique refuge. No Pleistocene refugia were previously identified in the central Andes for medium-sized mammals. However, considering the current distribution of the clades and the zones apparently less affected by the glaciations [5], clade A refuge probably existed in central Peru, while clade B and clade D refuges were possibly located in northern Chile and eastern Argentina, respectively. Clade C refuge was probably placed in eastern Bolivia. Individuals with haplotypes of clades B, C and D would subsequently migrate to a contact zone between 20 and 23°S during an interglacial period.

Implications for conservation

The concepts of "evolutionarily significant units" (ESUs) [41] and "management units" (MUs) [42] were created with the objective of targeting operational units for conservation below the species level. Although the use of the ESUs concept and its importance for prioritising conservation efforts is currently well accepted, its definition was strongly debated [43]. Many authors consider reciprocal monophyly as mandatory for the recognition of ESUs, as proposed by Moritz [42], or subspecies, while for other definitions a different evolutionary history between populations is sufficient [41,44-46]. MUs, in the other hand, are defined as population units with statistically significant differences of allelic frequencies at nuclear or mitochondrial level [42].

The long divergence time among clades, as well as the differentiation at the ncDNA level, highlight the importance of recognizing the four pampas cat population groups identified here (localities 2–10/11, 13–14/15–18/northern Chile) as different units for conservation in the Andean region. Since they have different allelic frequencies at both, mtDNA and ncDNA levels, all of these groups must to be recognised as MUs. In addition, we recommend the recognition of the localities 1, 12 and 19 as provisional MUs, since only private haplotypes were found for them (localities 1 and 19) or they differ at ncDNA level from adjacent MUs (locality 12; Figure 8).

The existence of an admixture zone in the central Andean region results in a lack of reciprocal monophyly between all the four population groups and makes their recognition as subspecies or ESUs controversial. However, cases of hybridization between subspecies or species are common in nature [47], of particular interest being the identification of hybridization areas for tigrinas, Geoffroy's cats and pampas cats [19,48]. In this sense, pampas cats from central Andes can be regarded as a complex of at least four ESUs/subspecies, one of them being the result of a past hybridization event.

Conclusion

Through the analyses of both mitochondrial and nuclear DNA, we inferred the population structure of the pampas cat in a broad portion of the Andean region, showing the existence of at least four groups of populations that must be considered different management units, and three other localities proposed here as provisional management units. The results revealed the influence of ancient climate fluctuations on the evolutionary history of this species, suggesting the split of four mtDNA clades during the Pre-Pastonian glacial period, long term isolation and further population expansions during the Aftonian interglacial. The pampas cat is a species with a relatively high mobility. It is expected that the other species of land vertebrates with a smaller or similar mobility have been affected in the same manner by the glacial and interglacial periods in the central Andes and, hence, they show similar population structures at present, that should be considered for the future study and conservation of the Andean fauna.

Methods

Study localities and sampling

A total of 39 skin samples and 532 faecal samples from 19 localities were analysed (Figure 1; Table 1). Samples were collected in the Andes of Argentina, Bolivia and Peru, covering more than 4350 km from 06°13'S to 41°07'S, with the exception of 5 tissue samples collected in the Argentine Pampas lowland of Buenos Aires province. Faecal samples were collected opportunistically in the field and kept in paper bags surrounded by silica gel in excess or in vials with ethanol for their transportation to the laboratory [49]. Once in the laboratory, samples were kept at -20°C, without silica gel, until DNA extraction. Skin samples were taken from stuffed animals owned by villagers or from museum specimens. For skin samples, approximately 0.5–1 cm2 was cut from the ear or, if that was not possible, from other areas, and kept in individual paper bags, in dry and cool conditions [50].

DNA extraction

DNA from skin samples was isolated by the standard method of proteinase K digestion, phenol-chloroform extraction and precipitation with ethanol [51]. DNA from faeces was isolated with the QIAamp® DNA Stool Mini Kit (QIAGEN, Ontario, Canada) according to the manufacturer's instructions, with the following two modifications: i) instead of 180–220 mg, from 200 to 500 mg of the superficial faecal material was used per sample, and ii) the ASL buffer was heated to 70°C before being added to the samples. To prevent and monitor the contamination of the samples during the laboratory processes, pre-PCR and post-PCR activities were carried out in different laboratories and negative controls were included in each batch of extraction and amplification [52,53].

Species identification

Faecal samples were identified to the species level by a PCR-RFLP method [54]. Briefly, a segment of 257–263 bp of the 16S mitochondrial gene was amplified by PCR and the product was exposed to the action of several restriction enzymes, resulting in species-specific fragment profiles which can be visualised on agarose gel.

MtDNA

The hypervariable domain 1 (HVS-I) of the mtDNA control region was selected because it is a non-coding sequence that is expected to display high polymorphism within felid species [27,55-57]. Because DNA from faeces and ancient tissues is often degraded [58,59] and prevents the amplification of the complete HVS-I, primers were designed to amplify two short segments on the HVS-I. In order to design these primers, the complete HVS-I region was amplified using fresh tissues from four Peruvian pampas cats with the primers CH3F and CH3R developed by Freeman et al. [56]. Both strands were sequenced with a CEQ 8000XL DNA Analysis System (Beckman Coulter Inc., Fullerton, Calif.). These sequences were aligned with the HVS-I region of domestic cat (Gene Bank accession number NC_001700), cheetah Acinonyx jubatus (NC_00512), margay Leopardus wiedii (AF129663S1-2) and ocelot Leopardus pardalis (AF129645S1-2) using the program Clustal W [60]. Finally, sequences conserved among species were used to design the felid-specific primers H1rev (5'-CCTGTACATGCTTAATATTC-3') and H2for (5'-ACATAYTATGTATATCGTGC-3') which provide PCR products smaller than 300 bp when used with CH3F and CH3R primers, respectively (Figure 2).

Amplification reactions were carried out in a volume of 12.5 μl containing a final concentration of 20 mM Tris-HCl (pH 8.4), 50 mM KCl, 1.5 mM MgCl2, 0.1 mM of each dNTP and 0.8 pM of each primer, 0.8 mg/ml of BSA, 0.2 unit of Taq DNA polymerase, and approximately 20 ng of template DNA. PCR conditions included an initial denaturing step at 92°C for 2 min, 45 cycles of 92°C for 15 s, 52°C for 15 s, and 68°C for 30 s, and a final extension step at 68°C for 5 min. Mitochondrial DNA variation was detected using single-strand conformation polymorphism (SSCP) [61]. The amplified products were electrophoresed on a 6% nondenaturing gel for 11 h 30' at 20 W in 0.5× TBE [62] and visualized using silver nitrate staining [63]. A band was sliced from the gel for each of the different observed conformers and placed in a volume of 35 μl of HPLC-grade water overnight. 2 μl of the dissolved DNA were used for amplification, and the products were sequenced. Because H1rev and H2for primers fit on a sequence of repeated tandems, called RS2 [28] (Figure 2), the product amplified with the CH3F-H1rev pair of primers was sequenced only in the forward direction, while that amplified with the H2for-CH3R primers was sequenced only in reverse. When possible, at least two individuals from different populations were sequenced for each observed conformer, to confirm the reliability of the SSCP protocol.

The mitochondrial genes NADH-5 (318 bp) and ATP-8 (191 bp) were sequenced for a sub-sample of 19 individuals, using the primers developed by Johnson et al. [29], to compare the phylogenetic relationships of the sampled individuals with pampas cats from other regions [20,29]. HVS-I, NADH-5 and ATP-8 sequences have been deposited in GeneBank (accession numbers FJ648644FJ648684, FJ664428FJ664446 and FJ664409FJ664427, respectively).

Microsatellites

Five microsatellite loci isolated from domestic cats (Fca24, Fca31, Fca45, Fca96 and Fca294) [64] were amplified and screened on acrylamide gels. Amplification reactions were carried out for all the pampas cat samples, with the same PCR conditions indicated above for the mtDNA. For the faecal samples, amplification and screening was made three times, following the multiple tube approach [65] and only the samples showing concordant results were used for data analysis.

Data analysis

Measures of population genetic diversity, including nucleotide diversity (π) and haplotype diversity (hd) [66] were estimated from mtDNA sequences using the computer program ARLEQUIN 2.00 [67]. Exact tests of population differentiation [68] with a 10000 Markov chain, and estimations of FST between pairs of localities [69] were performed with the same software.

The phylogenetic relationships of the aligned sequences were analysed by distance (neighbour joining, NJ) [70] and maximum likelihood (ML) methods using the program PAUP 4.01b [71]. Gaps were treated as single haplotype variants. After an election using MODELTEST 3.7 [72], the Tamura-Nei model [73] was used for NJ and ML analyses, assuming a gamma distribution for substitution rates across sites (AIC = 2534.70283). Node support was assessed using 1000 bootstrap replicates. Sequences of the equivalent control region segment from two ocelots and two margays were used as outgroups. To assess the extent of differentiation among regions and sampled localities, we performed an analysis of molecular variance (AMOVA) [74], with Φ-statistics, using ARLEQUIN 2.00 [67], with statistical significance tested using 10000 permutations. Several groupings were tested with AMOVA, and the grouping that maximized differentiation among regions and minimized differentiation among localities within regions was assumed to represent to the most parsimonious geographical subdivisions.

To estimate the divergence times between different groups, we considered the HVS-I region, a divergence time of 2.91 MY between the pampas cat and the ocelot and the margay, with a confidence interval between 2.02 and 4.25 MY [40], and the mean Tamura-Nei molecular distances between haplotypes of the identified clades as calculated in ARLEQUIN 2.00. All the identified haplotypes were used for each clade. Regions corresponding to deletions or insertions in different species were not considered for this analysis.

For each microsatellite genotype, the probability that two samples with the same genotype represent two different individuals was calculated from the allele frequencies in each sampled locality [75]. To obtain the value for the complete genotype, probabilities for each locus were multiplied, assuming independence of loci, as supported by the linkage map of microsatellite loci in the domestic cat [64].

Measures of expected heterozygosity (HE) and number of alleles (k) were estimated using the computer program ARLEQUIN 2.00 [67]. Hardy-Weinberg equilibrium was tested for each sampled locality with the method developed by Guo & Thompson [76] using GENEPOP 3.4 [77]. The population structure was examined for the microsatellite data by a Bayesian clustering method, using STRUCTURE 2.2 [78], with a 100000 burn-in period and 50000 MCMC repetitions, to infer the number of populations (K) and to assign individuals to inferred population clusters. Additionally, a neighbour-joining phylogenetic tree of the sampled localities with five or more individuals was constructed with the program POPULATIONS 1.2.28 [79] with the microsatellite data, using Dce distances [80] and 200 bootstraps on locus. As for the mtDNA data, AMOVA was performed on several groupings to test the population structure.

Authors' contributions

DC conceived the study and participated in its design, carried out the molecular genetic studies and the statistical analysis, and drafted the manuscript. BA coordinated the study, participated in the design and in the statistical analysis, and helped to draft the manuscript. ML and MRG contributed in data acquisition and interpretation, and helped to draft the manuscript. All authors read and approved the final manuscript.

Acknowledgements

Samples were generously provided by Susan Walker, Fernando Alfaro, Giovanna Gallardo, Rocío Palacios, Claudia Manfredi, Pablo Perovic, Lilian Villalba, Never Bonino, Victor Pacheco, Ursula Fajardo, Analí Madrid and José Luis Condori. We also would like to thank Mathieu Alday, Michelle Pelletier, Philippe Girard, Rachel Massicotte, Frederic Cyr, Claude-Olivier Silva-Beaudry, Émilie Castonguay, Joelle Boizard, Marie-Claire Binet and Juan Repucci, for their help during this project and to the Peruvian National Institute of Natural Resources – INRENA for research permits (N°13-2007) and institutional support. This work was supported by grants from Wildlife Conservation Network and the Panthera/Wildlife Conservation Society Kaplan Awards Program. Field collection of Argentine samples was supported by BP Conservation Programme, Whitley Fund for Nature, Wildlife Conservation Network and Darwin Initiative.

References

  1. Avise JC, Walker D, Johns GC: Speciation durations and Pleistocene effects on vertebrate phylogeography.

    P Roy Soc Lond B Bio 1998, 265:1707-1712. Publisher Full Text OpenURL

  2. Carstens BC, Brunsfeld SJ, Demboski JR, Good JD, Sullivan J: Investigating the evolutionary history of the Pacific Northwest mesic forest ecosystem: hypothesis testing within a comparative phylogeographic framework.

    Evolution 2005, 59:1639-1652. PubMed Abstract | Publisher Full Text OpenURL

  3. Hewitt GM: Some genetic consequences of ice ages, and their role in divergence and speciation.

    Biol J Linn Soc 1996, 58:247-276. OpenURL

  4. Hewitt GM: The genetic legacy of the Quaternary ice ages.

    Nature 2000, 405:907-913. PubMed Abstract | Publisher Full Text OpenURL

  5. Clapperton CM: Glacier Readvances in the Andes at 12500-10000 Yr Bp – Implications for Mechanism of Late-Glacial Climatic-Change.

    Journal of Quaternary Science 1993, 8(3):197-215. Publisher Full Text OpenURL

  6. Lavenu A: Formation and geological evolution. In Lake Titicaca: a synthesis of limnological knowledge. Edited by Dejoux C, Iltis A. Dordrecht, Nederlands: Kluwer Academic Publishers; 1992:3-15. OpenURL

  7. Markraf V: Climatic history of Central and South America since 18000 yr B.P: Comparison of pollen records and model simulations. In Global climates since the Last Glacial Maximum. Edited by Wrigth HE, Kutzbach JE, Webb T, Ruddiman WF, Street-Perrot FA, Bartlein PJ. Minnesota: University of Minnesota Press; 1993:357-385. OpenURL

  8. Nilsson T: The Pleistocene, geology and life in the quaternary ice age. Dordrecht, Nederlands: D. Reidel Publishing Company; 1983.

  9. Flenley JR: Tropical forests under the climates of the last 30,000 years.

    Climatic Change 1998, 39:177-197. Publisher Full Text OpenURL

  10. Graham A, Gregory-Wodzicki KM, Wright KL: Studies in Neotropical botany. XV. A Mio-Pliocene palynoflora from the Eastern Cordillera, Bolivia: implications for the uplift history of the Central Andes.

    Am J Bot 2001, 88:1545-1557. Publisher Full Text OpenURL

  11. Noonan BP, Gaucher P: Phylogegraphy and demography of Guianan harlequin toads (Atelopus), diversification within a refuge.

    Mol Ecol 2005, 14:3017-3031. PubMed Abstract | Publisher Full Text OpenURL

  12. Richardson JE, Pennington RT, Pennington TD, Hollingsworth PM: Rapid diversification of a species-rich genus of Neotropical rain forest trees.

    Science 2001, 293:2242-2245. PubMed Abstract | Publisher Full Text OpenURL

  13. Steele CA, Storfer A: Coalescent-based hypothesis testing supports multiple Pleistocene refugia in the Pacific Northwest for the Pacific giant salamander (Dicamptodon tenebrosus).

    Mol Ecol 2006, 15:2477-2487. PubMed Abstract | Publisher Full Text OpenURL

  14. Ursenbacher S, Carlsson M, Helfer V, Tegelstörm H, Fumagalli L: Phylogeography and Pleistocene refugia of the adder (Vipera berus) as inferred from mitochondrial DNA sequence data.

    Mol Ecol 2006, 15:3425-3437. PubMed Abstract | Publisher Full Text OpenURL

  15. Nowell K, Jackson P, IUCN/SSC Cat Specialist Group: Wild cats: status survey and conservation action plan. Gland, Switzerland: IUCN; 1996.

  16. Silveira L: Notes on the Distribution and Natural-History of the Pampas Cat, Felis-Colocolo, in Brazil.

    Mammalia 1995, 59(2):284-288. OpenURL

  17. García-Perea R: The pampas cat group (Genus Lynchailurus Severtzov, 1858) (Carnivora: Felidae), a systematic and biogeographic review.

    Am Mus Novit 1994, 3096:1-35. OpenURL

  18. Wozencraft WC: Order Carnivora. In Mammal species of the world. Edited by Wilson DE, Reeder DM. Baltimore and London: John Hopkins University Press; 2005:532-628. OpenURL

  19. Johnson WE, Slattery JP, Eizirik E, Kim JH, Raymond MM, Bonacic C, Cambre R, Crawshaw P, Nunes A, Seuanez HN, et al.: Disparate phylogeographic patterns of molecular genetic variation in four closely related South American small cat species.

    Molecular Ecology 1999, 8(12):S79-S94. PubMed Abstract | Publisher Full Text OpenURL

  20. Napolitano C, Bennett M, Johnson WE, O'Brien SJ, Marquet PA, Barria I, Poulin E, Iriarte A: Ecological and biogeographical inferences on two sympatric and enigmatic Andean cat species using genetic identification of faecal samples.

    Molecular Ecology 2008, 17(2):678-690. PubMed Abstract | Publisher Full Text OpenURL

  21. Nowell K: Revision of the Felidae Red List of Threatened Species.

    Cat News 2002, 37:4-6. OpenURL

  22. Cossíos ED, Madrid A, Condori JL, Fajardo U: An update on the distribution of Andean cat Oreailurus jacobita and pampas cat Lynchailurus colocolo in Peru.

    Endangered Species Research 2007, 3:313-320. Publisher Full Text OpenURL

  23. Villalba L, Lucherini M, Walker S, Cossíos D, Iriarte A, Sanderson J, Gallardo G, Alfaro F, Napolitano C, C S-Z: The Andean cat: A conservation action plan. La Paz, Bolivia: Andean Cat Alliance; 2004.

  24. Nowell K: The Cat Specialist Group digital library as a measure of cat conservation effort.

    Cat News 2002, 37:23-24. OpenURL

  25. de Oliveira T, Eizirik E, Lucherini M, Acosta G, Leite-Pitman R, Pereira J: Leopardus colocolo. [http://www.iucnredlist.org] webcite

    In 2008 IUCN Red List of Threatened Species Edited by IUCN. 2008.

    Downloaded on 14 January 2009

    OpenURL

  26. Burger PA, Steinborn R, Walzer C, Petit T, Mueller M, Schwarzenberger F: Analysis of the mitochondrial genome of cheetahs (Acinonyx jubatus) with neurodegenerative disease.

    Gene 2004, 338(1):111-119. PubMed Abstract | Publisher Full Text OpenURL

  27. Eizirik E, Bonatto SL, Johnson WE, Crawshaw PG, Vie JC, Brousset DM, O'Brien SJ, Salzano FM: Phylogeographic patterns and evolution of the mitochondrial DNA control region in two neotropical cats (Mammalia, Felidae).

    Journal of Molecular Evolution 1998, 47(5):613-624. PubMed Abstract | Publisher Full Text OpenURL

  28. Lopez JV, Cevario S, OBrien SJ: Complete nucleotide sequences of the domestic cat (Felis catus) mitochondrial genome and a transposed mtDNA tandem repeat (Numt) in the nuclear genome.

    Genomics 1996, 33(2):229-246. PubMed Abstract | Publisher Full Text OpenURL

  29. Johnson WE, Culver M, Iriarte JA, Eizirik E, Seymour KL, O'Brien SJ: Tracking the evolution of the elusive Andean mountain cat (Oreailurus jacobita) from mitochondrial DNA.

    Journal of Heredity 1998, 89(3):227-232. PubMed Abstract | Publisher Full Text OpenURL

  30. Goudet J: FSTAT (vers. 1.2): a computer program to calculate F-statistics.

    J Hered 1995, 86:485-486. OpenURL

  31. Bamshad MJ, Wooding S, Watkins WS, Ostler CT, Batzer MA, Jorde LB: Human population genetic structure and inference of group membership.

    American Journal of Human Genetics 2003, 72:578-589. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  32. Ammann C, Jenny B, Kammer K, Messerli B: Late Quaternary Glacier response to humidity changes in the arid Andes of Chile (18–29 degrees S).

    Palaeogeography Palaeoclimatology Palaeoecology 2001, 172(3–4):313-326. Publisher Full Text OpenURL

  33. Kull C, Grosjean M, Veit H: Modeling modern and Late Pleistocene glacio-climatological conditions in the north Chilean Andes (29–30 degrees S).

    Climatic Change 2002, 52(3):359-381. Publisher Full Text OpenURL

  34. Marin JC, Casey CS, Kadwell M, Yaya K, Hoces D, Olazabal J, Rosadio R, Rodriguez J, Spotorno A, Bruford MW, et al.: Mitochondrial phylogeography and demographic history of the Vicuna: implications for conservation.

    Heredity 2007, 99(1):70-80. PubMed Abstract | Publisher Full Text OpenURL

  35. Yensen E, Tarifa T: Galictis cuja.

    Mam Spec 2003, 728:1-8. Publisher Full Text OpenURL

  36. Sheffield SR, Thomas HH: Mustela frenata.

    Mammalian Species 1997, 570:1-9. Publisher Full Text OpenURL

  37. Ruiz-García M, Pacheco L, Alvarez D: Caracterización genética del puma andino boliviano.

    Revista Chilena de Historia Natural 2008, 81(4):470-492. OpenURL

  38. Mercer JH: Glacial history of Southernmost South America.

    Quaternary Res 1976, 6:125-166. Publisher Full Text OpenURL

  39. Chesser RT: Evolution in the high Andes: the phylogenetics of Muscisaxicola groud-tyrants.

    Mol Phylogenet Evol 2000, 15:369-380. PubMed Abstract | Publisher Full Text OpenURL

  40. Johnson WE, Eizirik E, Pecon-Slattery J, Murphy WJ, Antunes A, Teeling E, O'Brien SJ: The Late Miocene radiation of modern Felidae: A genetic assessment.

    Science 2006, 311(5757):73-77. PubMed Abstract | Publisher Full Text OpenURL

  41. Ryder OA: Species Conservation and Systematics – the Dilemma of Subspecies.

    Trends in Ecology & Evolution 1986, 1(1):9-10. Publisher Full Text OpenURL

  42. Moritz C: Defining Evolutionarily-Significant-Units for Conservation.

    Trends in Ecology & Evolution 1994, 9(10):373-375. Publisher Full Text OpenURL

  43. Fraser DJ, Bernatchez L: Adaptive evolutionary conservation: towards a unified concept for defining conservation units.

    Molecular Ecology 2001, 10(12):2741-2752. PubMed Abstract | Publisher Full Text OpenURL

  44. Dizon AE, Lockyer C, Perrin WF, Demaster DP, Sisson J: Rethinking the Stock Concept – a Phylogeographic Approach.

    Conservation Biology 1992, 6(1):24-36. Publisher Full Text OpenURL

  45. Vogler AP, DeSalle R: Diagnosing units of conservation management.

    Conserv Biol 1994, 6:170-178. OpenURL

  46. Waples RS: Pacific Salmon, Oncorhynchus spp. & the definition of 'species' under the endangered species act.

    Mar Fish Rev 1991, 53:11-22. OpenURL

  47. Arnold ML: Natural hybridization and evolution. Oxford: Oxford University Press; 1997.

  48. Eizirik E, Indrusiak C, Trigo TC, Sana DA, Mazim FD, Freitas TRO: Refined Mapping and Characterization of the Geographic Contact Zone Between Oncilla and Geoffroy's Cat in Southern Brazil.

    IUCN Cat News 2006, 45:8-11. OpenURL

  49. Wasser SK, Houston CS, Koehler GM, Cadd GG, Fain SR: Techniques for application of faecal DNA methods to field studies of Ursids.

    Molecular Ecology 1997, 6(11):1091-1097. PubMed Abstract | Publisher Full Text OpenURL

  50. Cossíos D, Beltrán Saavedra F, Bennet M, Bernal N, Fajardo U, Lucherini M, Merino J, Marino J, Napolitano C, Palacios R, et al.: Manual de metodologías para relevamientos de carnívoros altoandinos. Buenos Aires, Argentina: Alianza Gato Andino; 2007.

  51. Sambrook J, Fritsch EF, Maniatis T: Molecular cloning: a laboratory manual. 2nd edition. Cold Spring Harbor, N.Y.: Cold Spring Harbor Laboratory; 1989.

  52. Kohn MH, Wayne RK: Facts from feces revisited.

    Trends in Ecology & Evolution 1997, 12(6):223-227. Publisher Full Text OpenURL

  53. Taberlet P, Waits LP, Luikart G: Noninvasive genetic sampling: look before you leap.

    Trends Ecol Evol 1999, 14:323-327. PubMed Abstract | Publisher Full Text OpenURL

  54. Cossíos ED, Angers B: Identification of Andean felid feces using PCR-RFLP.

    Journal of Neotropical Mammalogy 2006, 13:239-244. OpenURL

  55. Eizirik E, Kim JH, Menotti-Raymond M, Crawshaw PG, O'Brien SJ, Johnson WE: Phylogeography, population history and conservation genetics of jaguars (Panthera onca, Mammalia, Felidae).

    Molecular Ecology 2001, 10(1):65-79. PubMed Abstract | Publisher Full Text OpenURL

  56. Freeman AR, MacHugh DE, McKeown S, Walzer C, McConnell DJ, Bradley DG: Sequence variation in the mitochondrial DNA control region of wild African cheetahs (Acinonyx jubatus).

    Heredity 2001, 86:355-362. PubMed Abstract | Publisher Full Text OpenURL

  57. Johnson WE, Godoy JA, Palomares F, Delibes M, Fernandes M, Revilla E, O'Brien SJ: Phylogenetic and phylogeographic analysis of Iberian lynx populations.

    Journal of Heredity 2004, 95(1):19-28. PubMed Abstract | Publisher Full Text OpenURL

  58. Frantzen MAJ, Silk JB, Ferguson JWH, Wayne RK, Kohn MH: Empirical evaluation of preservation methods for faecal DNA.

    Molecular Ecology 1998, 7(10):1423-1428. PubMed Abstract | Publisher Full Text OpenURL

  59. Lindahl T: Instability and Decay of the Primary Structure of DNA.

    Nature 1993, 362(6422):709-715. PubMed Abstract | Publisher Full Text OpenURL

  60. Thompson JD, Higgins DG, Gibson TJ: CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice.

    Nucleic Acids Res 1994, 22:4673-4680. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  61. Orita M, Iwahana H, Kanazawa H, Hayashi K, Sekiya T: Detection of Polymorphisms of Human DNA by Gel-Electrophoresis as Single-Strand Conformation Polymorphisms.

    Proceedings of the National Academy of Sciences of the United States of America 1989, 86(8):2766-2770. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  62. Angers B, Bernatchez L: Combined use of SMM and non-SMM methods to infer fine structure and evolutionary history of closely related brook charr (Salvelinus fontinalis, Salmonidae) populations from microsatellites.

    Molecular Biology and Evolution 1998, 15(2):143-159. OpenURL

  63. Bassam BJ, Caetanoanolles G, Gresshoff PM: Fast and Sensitive Silver Staining of DNA in Polyacrylamide Gels.

    Analytical Biochemistry 1991, 196(1):80-83. PubMed Abstract | Publisher Full Text OpenURL

  64. Menotti-Raymond M, David VA, Lyons LA, Schaffer AA, Tomlin JF, Hutton MK, O'Brien SJ: A genetic linkage map of microsatellites in the domestic cat (Felis catus).

    Genomics 1999, 57(1):9-23. PubMed Abstract | Publisher Full Text OpenURL

  65. Taberlet P, Griffin S, Goossens B, Questiau S, Manceau V, Escaravage N, Waits LP, Bouvet J: Reliable genotyping of samples with very low DNA quantities using PCR.

    Nucleic Acids Research 1996, 24(16):3189-3194. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  66. Nei M: Molecular evolutionary genetics. New York: Columbia University Press; 1987.

  67. Schneider S, Roessli D, Excofier L: Arlequin: A software for population genetics data analysis. 2.000th edition. Geneva: Genetics and Biometry Lab, Dept. of Anthropology, University of Geneva; 2000.

  68. Raymond M, Rousset F: An exact test for population differentiation.

    Evolution 1995, 49(6):1280-1283. Publisher Full Text OpenURL

  69. Weir BS, Cockerham CC: Estimating F-statistics for the analysis of population structure.

    Evolution 1984, 38:1358-1370. Publisher Full Text OpenURL

  70. Saitou N, Nei M: The Neighbor-Joining Method – a New Method for Reconstructing Phylogenetic Trees.

    Mol Biol Evol 1987, 4(4):406-425. PubMed Abstract | Publisher Full Text OpenURL

  71. Swofford DL: PAUP*: phylogenetic analysis using parsimony (and other methods). 4.06b edition. Massachusetts: Sinauer Associates; 1993. PubMed Abstract OpenURL

  72. Posada D, Crandall KA: MODELTEST: testing the model of DNA substitution.

    Bioinformatics 1998, 14(9):817-818. PubMed Abstract | Publisher Full Text OpenURL

  73. Tamura K, Nei M: Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees.

    Mol Biol Evol 1993, 10:512-526. PubMed Abstract | Publisher Full Text OpenURL

  74. Excoffier L, Smouse PE, Quattro JM: Analysis of Molecular Variance Inferred from Metric Distances among DNA Haplotypes – Application to Human Mitochondrial-DNA Restriction Data.

    Genetics 1992, 131(2):479-491. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  75. Paetkau D, Strobeck C: Microsatellite Analysis of Genetic-Variation in Black Bear Populations.

    Molecular Ecology 1994, 3(5):489-495. PubMed Abstract | Publisher Full Text OpenURL

  76. Guo SW, Thompson EA: Performing the Exact Test of Hardy-Weinberg Proportion for Multiple Alleles.

    Biometrics 1992, 48(2):361-372. PubMed Abstract | Publisher Full Text OpenURL

  77. Raymond M, Rousset F: Genepop (Version-1.2) – Population-Genetics Software for Exact Tests and Ecumenicism.

    Journal of Heredity 1995, 86(3):248-249. OpenURL

  78. Pritchard JK, Stephens M, Donnelly P: Inference of population structure using multilocus genotype data.

    Genetics 2000, 155(2):945-959. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  79. Langella O: POPULATIONS, a free population genetic software.

    1.2.28 (12/5/2002) edn: CNRS UPR9034 2002. OpenURL

  80. Cavalli-Sforza LL, Edwards AWF: Phylogenetic analysis: models and estimation procedures.

    Evolution 1967, 32:550-570. Publisher Full Text OpenURL