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

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][2][3][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.

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.
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][27][28].
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.
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.

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.
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)(16)(17)(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).
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 Total

(39) 199
After the total number of pampas cat samples, number of skin samples is indicated in parentheses. Number of individuals as identified from microsatellite data.
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.

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).
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.
Pampas cat sampled localities Figure 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).
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.  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  [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], H1rev H2for CHR3 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 con-servation 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][45][46]. MUs, in the other hand, are defined as population units with statistically significant differences of allelic frequencies at nuclear or mitochondrial level [42].
Neighbour-joining tree of the observed HVS-Isequences Figure 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.
Neighbour-joining tree of the combined observed NADH-5 -ATP-8 sequences 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 pro-visional 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 com- The locality numbers correspond to those on the Table 1. Mean values were calculated with the localities with more than 4 individuals only. Total values were calculated for the entire sample. n = sample size; H = number of haplotypes; hd = haplotype diversity; S = number of polymorphic sites; P = percentage of variable sites; π = nucleotid diversity; k = number of alleles; H E = expected heterozygosity. Analyses were performed on mtDNA haplotypes using Φ-statistics or microsatellite data using R ST -statistics and distances. Localities with less than 5 individuals were not included. mon 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.

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, cov-Population structure inferred from microsatellite data   ering 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, approxi-mately 0.5-1 cm 2 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][56][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 (   [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.

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 F ST 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 (H E ) 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 STRUC-TURE 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