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 FJ648644–FJ648684, FJ664428–FJ664446 and FJ664409–FJ664427, 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.