Sampling
Scat collection was conducted by Panthera Belize in two sites (Figure 2): 1) the Cockscomb Basin Wildlife Sanctuary and some areas that fall beyond the boundaries of this natural reserve, between 2003-2007; and 2) the Maya Forest Corridor between 2009-2011. Sampling was conducted as part of the Global Felid Genetics Program. The Maya Forest Corridor (17.349140° N, 88.455310° W, 50m elevation) has been identified as an important link between jaguar populations in northern and southern Belize (20). The Cockscomb Basin Wildlife Sanctuary (16.7162° N, 88.6608° W, 500 m elevation) supports the highest density of jaguars in Mesoamerica (15,19). These two areas play a crucial role in the maintenance of the Mesoamerican Biological Corridor, comprised of a network of protected areas stretching from Mexico to Panama (23). A total of 852 scats were collected opportunistically along forest trails and unpaved roads. Samples were georeferenced at the time of collection and stored individually at room temperature using silica gel beads until DNA extractions. All faecal material has been deposited in the Sackler Institute for Comparative Genomics at the American Museum of Natural History, NYC.
DNA extraction and species identification
About 200mg of the dry sample was shaved from each scat and used to extract genomic DNA using a QIAmp DNA extraction Stool Mini Kit (Qiagen Inc., Valencia, California, USA). Samples were screened to identify species via PCR using three mitochondrial gene regions including cytochrome b (H15149)(49,50), 12S rDNA (L1085 and H1259)(51) 16S rDNA (L2513 and H2714)(51) and 16Scp (16S cp-F 16S cp-R)(51). PCR amplifications were done using G&E Ready-to-go PCR Beads (GE Healthcare, Piscataway, New Jersey, USA) in 25ul reaction volumes containing 0.2uM of each primer, 0.3uL BSA and 1ul of DNA. PCR profiles for each reaction are available as supplementary material (see Additional file 2). All extractions and amplifications included negative controls. Cycle-sequencing was performed with the BidDye® Terminator v. 1.1 cycle Sequencing kit (Applied Biosystems, Lennik, Belgium) using the PCR primers. Amplified products were visualised and scored in an ABI 3730xl DNA Analyzer (Applied Biosystems, Carlsbad, CA). Sequences were aligned and manually corrected using GENEIOUS v.6.5 (Biomatters Ltd., Auckland, New Zealand). Species identification was confirmed by comparing consensus sequences to known felid references and by constructing a phylogenetic tree with 93% similarity using the Jukes-Cantor neighbourhood-joining model (52). Species ID was validated if at least three gene fragments successfully identified the same species.
Individual and sex identification
Samples identified as P. onca were genotyped for twelve polymorphic microsatellite loci (FCA 32, 75, 96, 126, 100, 124, 132, 208, 212, 229, 275, and 225) (53). Each sample was screened 4 times to have enough replicates to assign the alleles correctly and to reduce artefacts caused by allele dropout, very difficult samples were screened up to 6 times. The genotypes were validated if they successfully scored at least three times. Samples that failed to amplify at least six loci were not considered in our analysis. We amplified sets of 2-3 primers in multiplex PCR reactions and FCA225 as singleplex (see Additional file 2). Amplification was performed using the Qiagen Multiplex Master Mix ® (Qiagen Inc., Valencia, CA), with a final volume of 20ul following manufacturers recommendations.
PCR reaction conditions were optimised in the annealing temperature for each set of primers as follows: 95°C for 15min, followed by 13 cycles of 94 °C for 30s, 57.4-62.4 °C for 90s, 72°C for 60s, followed by 32 cycles of 94°C for 30s, 55-60 °C for 90 s, 72°C for 60s and 60 °C for 30min (see Additional file 2). Amplified products were visualised and scored in an ABI 3730xl DNA Analyzer (Applied Biosystems, Carlsbad, CA). PCR products were analysed using GeneMapper v5.0 following the quality index measure to validate genotypes (Miquel et al. 2006). In order to further validate the allele calling, we used FRAGMAN v.1.0.8 package in R (54). Identification of unique multilocus genotypes was achieved using the ALLELEMATCH package in R v.2.5 (55), which accommodates for genotyping error and missing data by making pairwise comparisons and finding similarity scores for each pair of profiles and clustering similar ones. To determine if genotypes were not different individuals matched by chance, the cumulative probability of identity for unrelated individuals (PID) and siblings (Psib) was calculated. Genotypes with Psib < 0.010 were positively identified as single individuals (25,56).
Gender was determined following the PCR-CTPP method for sex identification developed by Wei et al. (2008) based on zinc finger alleles (ZFX/ZFY). The two forward primers (ZF-1F and ZFY-2F) were fluorescently labelled with FAM and MAX, respectively. Amplification was performed using G&E Ready-to-go PCR Beads (GE Healthcare, Piscataway, New Jersey, USA) in a final volume of 25ul containing 2mM of each primer and 5ul of DNA. The PCR conditions were as follows: 95°C for 5min, followed by 13 cycles of 95 °C for 45s, 63 °C for 60s, 72°C for 60s, followed by 35 cycles of 95°C for 45s, 53 °C for 45s, 72°C for 60s and 72 °C for 10min (see Additional file 2). DNA amplifications were visualised in a 5% agarose gel to confirm the molecular sexing. Amplified products were visualised and scored in an ABI 3730xl DNA Analyzer (Applied Biosystems, Carlsbad, CA). PCR products were analysed using GeneMapper v5.0.
Conservation genetics analyses
MICROCHECKER v 2.2.3 (57) was used to screen null alleles at each locus. Measures of genetic diversity as average number of alleles per locus (A), allelic richness (AR), observed (Ho) and expected heterozygosity (He) deviations from Hardy-Weinberg equilibrium (HWE), inbreeding coefficient (FIS), and multivariate analyses were performed using the ADEGENET package in R v.2.0.1 (58). The rarefaction procedure implemented the HIERFSTAT package in R v 2.9.3.250 (59) was used to estimate the expected number of alleles (NA) and to compare allelic richness (AR). Genetic differentiation among sampled individuals was summarised in a Discriminant analysis of principal components (DAPC) based on the allele frequencies. Linkage disequilibrium (LD) between pairs of loci was performed in GENEPOP v. 4.2 (60) with default settings. We used GenAlex v6.0 (61) to test the occurrence of a positive correlation (Rxy>0) between the genetic PHIPT matrix and geographic distance so-called Isolation by Distance (IBD) via Mantel tests of matrix correspondence; we assessed the partitioning of genetic variation between sampling localities with Analysis of Molecular Variance (AMOVA).
Population Structure
We estimated population genetic structure using Bayesian assignment methods with STRUCTURE v2.3.4 (62), which assigns individuals to a number K of genetically homogeneous groups, based on the Bayesian clustering in accordance to the expected Hardy–Weinberg equilibrium and absence of linkage disequilibrium between loci. We ran STRUCTURE using the admixture model with correlated allele frequencies and sampling locations as spatial prior (LOCPRIOR) to allow sampling location to assist in the clustering, and we performed 20 independent runs for K = 1-10. We set a burn-in period of 100,000 and 1,000,000 MCMC iterations and assumed an admixture model with correlated allele frequencies. To determine the optimal number of clusters and render bar plots, we implemented the Evanno method (35) using POPHELPER package in R v1.2.1 (63). Furthermore, we inferred spatial genetic structure with TESS v2.3.1 (64). This program assumes that population memberships follow a hidden Markov random field model where the log-probability of an individual belonging to a particular population, given the population membership of its closest neighbours, is equal to the number of neighbours belonging to this population (65). We tested the CAR, and BYM models with linear trend surface to define the spatial prior for admixture (36); we set a burn-in period of 100,000 and 1000,000 sweeps through 10 independent runs testing the maximal number of clusters from 1 – 10. To decide the optimal K, we plotted the deviance information criterion (DIC) against Kmax. We also used a spatial clustering approach in GENELAND v4.0.6 (65) as an additional method to infer the number of populations and the spatial location of genetic discontinuities. This program allows using georeferenced individual multilocus genotypes to infer the number of populations and uses the spatial location of genetic discontinuities between those populations. We determined K across 20 independent runs with 1,000,000 MCMC iterations. Thinning was set at 100, allowing K to vary from 1 to 10. We used the correlated allele model and set the maximum rate of the Poisson process at 50 (the number of individuals), the maximum number of nuclei in the Poisson-Voronoi tessellation at 150 (three times the number of individuals), and the uncertainty of spatial coordinates of the collection at 25 meters. We re-ran the analysis ten times to check for consistency across runs.
To further explore genetic diversity and structure among individuals, we reduced the dimensions via a Discriminant Principal Component Analysis (DAPC) without a priori group assignment using the ADEGENET package in R v2.0.1 (58,66). The tools implemented in DAPC allow solving complex population structures by summarising the genetic differentiation between groups while overlooking within-group variation, therefore achieving the best discrimination of individuals into pre-defined groups (67). This multivariate method is useful to identify clusters of genetically related individuals when group priors are lacking. Estimation of clusters was performed by comparing the different clustering solutions using the Bayesian Information Criterion (BIC). We compared the results from the three Bayesian approaches and the DAPC to provide confidence in the spatial designation of genetic groupings.
Relatedness
Levels of genetic relatedness were calculated using seven estimators as implemented in the RELATED package in R v1.0 (68). Pairwise relatedness was calculated using the estimators described by Queller and Goodnight, 1989; Li et al., 1993; Ritland, 1996; Lynch and Ritland, 1999 and Wang, 2002), as well as the dyadic likelihood estimator described in Milligan (2003) and the triadic likelihood estimator from Wang (2007). Genotyping errors and inbreeding estimations were incorporated into the model, and confidence intervals (95%) were obtained through bootstrapping across loci. Allele frequencies were used to simulate pairs of individuals of known relatedness based on Parent-Offspring, Full siblings, Half siblings and Unrelated individuals.
Landscape Connectivity
We used ArcGIS Desktop v10 to model resistance surfaces of the landscape to jaguar movement and CIRCUITSCAPE v3.5 (76) to identify the most probable routes for dispersal and gene flow between the Cockscomb Basin Wildlife Sanctuary and the Maya Forest Corridor by. Our analysis describes the connectivity and permeability of the landscape between the Cockscomb Basin Wildlife Sanctuary and the Maya Forest and was used to predict all possible paths connecting these two areas. We created a resistance raster surface at a 50 metres resolution by combining data on ecosystem type, human settlement and roads. Spatial data were obtained from the Biodiversity and Environmental Resource Data System of Belize (BERDS) and included ecosystem types (2017), protected areas (2015), roads (2013), and human settlements (2014) (77).
The probability of jaguar occurrence has been positively associated with forest cover and negatively associated with human activities (78). We assumed that landscape features vary in their importance to promote connectivity for jaguars; for example, lowland broad-leaved forests may promote jaguar movement while human settlement may limit it. To represent the variation in importance for connecting different areas, we created a surface of cost values of habitat preference (Table 2) based on published scientific literature and expert knowledge of nine jaguar scientists with experience working in Belize and South-Eastern Mexico. Although basing cost values on expert knowledge as quantitative information is controversial, it has been commonly used as a surrogate when empirical data is limited or unavailable (8,75–78). To gather the information, we asked jaguar scientists to rank habitat preference according to their empirical and theoretical knowledge on habitat preference and resistance features to allow movement of jaguars across the landscape. The expert’s knowledge was supplemented with empirical data of jaguar habitat use (4,16,89,23,39,83–88). Therefore, by expert knowledge, we mean the experts’ personal judgement that agrees with empirical data.
Each pixel was assigned a cost value in a scale of 1 – 9, which represents the relative effort required to move from one point to another (higher values facilitate movement) (Table 2). The average locations of the identified jaguars were grouped in single polygons per study area and were calculated as the 95% minimum convex polygon (MCP). The resulting connectivity between polygons is represented as a cumulative current flow map where areas with higher predicted movement are coloured in red and areas of lower predicted movement in blue. The cumulative current maps represent all potential routes that an animal moving from one polygon could use to move across central Belize and potentially disperse to and from the Cockscomb Basin Wildlife Sanctuary to the Maya Forest Corridor (Figure 2).
Table 2. Cost-values of jaguar habitat preference based on expert knowledge. Ecosystem names follow UNESCO’s classification system. Values range from 1 (not preferred) to 9 (highly preferred).
|
|
Landscape feature
|
Preference
|
Lowland broad-leaved wet forest
|
9
|
Submontane broad-leaved wet forest
|
8.8
|
Lowland broad-leaved moist forest
|
8.7
|
Submontane broad-leaved moist forest
|
8.7
|
Lowland broad-leaved moist scrub forest
|
7.4
|
Lowland broad-leaved wet forest: Steep
|
7.1
|
Submontane broad-leaved wet forest: Steep
|
7
|
Lowland broad-leaved dry forest
|
7
|
Submontane broad-leaved moist forest: Steep
|
6.9
|
Lowland broad-leaved moist forest: Steep
|
6.9
|
Lowland broad-leaved moist scrub forest: Steep
|
6.1
|
Mangrove and littoral forest
|
5.8
|
Lowland pine forest
|
5.7
|
Wetland
|
5.4
|
Submontane pine forest
|
5.4
|
Water
|
5.3
|
Shrubland
|
5.1
|
Lowland savanna
|
4.4
|
Submontane pine forest: Steep
|
3.9
|
Anthropogenic feature
|
>5 km distance from Highways
|
9
|
>5 km distance from human settlements
|
9
|
> 3 and ≤ 5km distance from Highways
|
6.5
|
> 3 and ≤ 5km distance from human settlements
|
6.5
|
≤ 3km distance from Highways
|
6.1
|
≤ 3km distance from human settlements
|
5.4
|
Highway
|
3.3
|
Agricultural uses
|
2.7
|
Aquaculture
|
2
|
Wasteland
|
1.8
|
Urban
|
1
|