Population Genetic Structure and Habitat Connectivity for Jaguar (Panthera onca) Conservation in Belize


 Background: Effective connectivity between jaguar ( Panthera onca ) populations across the American continent will ensure the natural gene flow and the long-term survival of the species throughout its range. Jaguar conservation efforts have focused primarily on connecting suitable habitat in a broad-scale. However, accelerated habitat reduction, limited funding, and the complexity of jaguar behaviour have proven challenging to maintain connectivity between populations effectively. Here we used individual-based genetic analysis in synthesis with landscape permeability models to assess levels of current genetic connectivity and identify alternative corridors for jaguar movement between two core areas in central and southern Belize. Results: We use 12 highly polymorphic microsatellite loci to identify 50 distinct individual jaguars, including 41 males, 3 females and 6 undetermined animals, from scat samples collected in The Cockscomb Basin Wildlife Sanctuary and The Central Belize Corridor. Using Bayesian and multivariate models of genetic structure, we identified one single group across two sampling sites with low genetic differentiation between them. We used fine-scale data on biodiversity features as vegetation types to predict the most probable corridors using least-cost path analysis and circuit theory. Conclusions: The results of our study highlight the importance of expanding the boundaries of the Central Belize Corridor to effectively cover areas that would more easily facilitate jaguar movement between locations.

by correlating the spatial heterogeneity of landscapes with estimates of gene flow (18).
This approach provides important opportunities to identify areas of conservation priority and provide critical knowledge on habitat fragmentation, dispersal ecology, and functional connectivity in complex landscapes for effective conservation planning and management actions (19).
The establishment of corridors to improve population connectivity, particularly in Mesoamerica, has been among some of the most important efforts to prevent the loss of biodiversity in the world's biologically richest regions (20). One of the most important jaguar populations in Mesoamerica can be found in the forests of Belize, a core and critical area for the species throughout its range (21). According to the Belize National Protected Areas System Plan (22), 36% of Belize land territory is under conservation management. In particular, with its 390 square kilometres of protected forests, The Cockscomb Basin Wildlife Sanctuary harbours one of the largest concentrations of jaguars in the country (21,23). This area is connected to other Natural Protected Areas via the Central Belize Corridor, which extends over 750 square kilometres and is considered the most critical and important corridor of the Belize National Protected Areas Systems (24). Range-wide corridors have been established as a major tool to improve population connectivity and thus aid the subsistence of jaguars (5). However, to build a well-connected protection network, considerations must be taken on the spatial scale at which conservation strategies are implemented.
Broad-scale conservation efforts will benefit from information gained at a finer scale, especially in heterogeneous or fragmented areas (3, 9,25). The effective collaboration among scientists, practitioners, non-governmental organisations and politicians will tap the full potential of reservoir projects and conservation actions across the jaguar's range.
Here, we present a comprehensive study on genetic structuring and patterns of landscape 5 heterogeneity to identify alternative habitat corridors for gene flow between populations within two principal locations within a jaguar stronghold in Belize. Using jaguar scat samples, we investigate population genetic structure, levels of inbreeding and gene flow.
Additionally, we correlate landscape features and patterns of gene flow to examine landscape permeability for jaguars between The Cockscomb Basin Wildlife Sanctuary and The Central Belize Corridor.

Genetic variation
A total of 536 scat samples collected across the two sampling areas were positively matched to P. onca. Other identified species included Puma concolor, Leopardus wideii, Leopardus pardalis, Herpailurus yaguarondi, and Canis familiaris. Genotyping revealed a total of 50 unique multilocus genotypes (37 from the Cockscomb reserve and 13 from the Central Corridor); these included 41 prospective males, 3 prospective females and 6 unidentified genders (Additional file 3). MICROCHEKER detected three loci (FCA 212, 229, and 075), showing signs of a null allele, but did not find evidence of scoring mistakes or large allele dropout (Additional file 1).
The mean polymorphic information content PIC = 0.642. The genotyping results allowed the identification of 50 distinct individual jaguars, 37 corresponding to Cockscomb Reserve and 13 to Central Corridor, the geographical coordinates assigned to each individual was determined by averaging the coordinates of all the samples corresponding to that particular individual. Tests for departure from Hardy-Weinberg Equilibrium were 6 variable for each locus, with four loci deviating from HWE (Table1). This deviation could be explained by a deficit of heterozygotes within the population potentially caused by inbreeding or by the presence of null alleles (FIS = 0.22, p-value = 0.001). Furthermore, this test showed evidence of low genetic differentiation (FST = 0.021, p-value = 0.007; FIT = 0.237, p-value = 0.001); linkage disequilibrium was not significant for any pair of loci.
The PCA of genetic diversity shows overlapping of the two sites, indicating overlapping of allele frequencies and little differentiation between groups ( Figure 1D). The AMOVA analysis revealed that less than 2% of genetic variance occurred among individuals in the Cockscomb Reserve and individuals in the Central Corridor; and showed low levels of genetic differentiation between groups (FST = 0.015, p-value = 0.026). Results from the Mantel test showed significant evidence of isolation by distance (Rxy = 0.167, p-value = 0.01; Figure 1C).

Population Structure and relatedness
Data analysis using STRUCTURE revealed that k = 1 had the highest mean probability of density value, and k = 4 had the highest delta-K value ( Figure 1B). This was consistent with the results from TESS where k = 1 also had the highest probability (ΔK = 10.15), in both cases no clear pattern of genetic structure can be observed when rendering the assignment probability in bar plots ( Figure 1B). Results from GENELAND revealed that k = 1 alsohad the highest probability in 8 of 10 runs, and the final map does not show a clear population boundary between sampling sites ( Figure 1E). The DAPC analysis showed that the lowest BIC value (68.42) corresponded to K = 2 and is represented in a single discriminant function; however, the BIC difference between K = 2 and K = 1 is negligible The performance of the seven relatedness estimators was analysed to provide information on the degree of resolution expected in our dataset. Mean relatedness amongst 7 individuals from the Central Corridor was -0.046 ± 0.068 (SE = 0.008) and from Cockscomb Reserve -0.15 ± 0.086 (SE = 0.003). Amongst all individuals mean relatedness was -0.01 ± SD 0.08 (SE = 0.002). Overall, individuals from the Central Corridor were more closely related to each other than to those in the Cockscomb Reserve.

Landscape permeability
The least-cost path analysis inferred the best route from the most northern point in the

Population genetics
This study presents estimates of genetic variation for individuals sampled within two core areas in Central and Southern Belize. Twelve polymorphic microsatellite loci were useful to successfully identify 50 distinct individuals that correspond to 41 males, 3 females and 6 undetermined sexes. The relatively low number of females could be explained by the sampling methodology rather than reflecting the proportion of sexes in the area. Sampling was conducted close to paths or dirt roads, and closer to human settlements; because 8 females could be more elusive, have smaller home ranges, hide their scats and avoid crossing open spaces and wide paths (3,12,26) this method could favour the sampling of male scats and therefore bias the analysis towards the more frequently observed sex.
Studies on dispersal in large felines show that males are the dispersing sex, while females tend to be more philopatric (9,11,27,28); other measurements for genetic differentiation between sexes and more female scat samples are necessary to confirm sex-biased dispersal in this area.  (29). Additionally, properties of our data as sample size, number of loci, polymorphisms, and null alleles could have also influenced their performance (30,31). Their power to detect population structure has been shown to decrease in accuracy at very low levels of population differentiation (FST<0.05) (29) as the number of estimated populations can be affected by a violation of model assumptions and cryptic relatedness (32). Furthermore, GENELAND and DAPC also indicated a single cluster as the most probable number of populations. GENELAND assigned individuals to each population considering the sampling locations and measurements of genetic 9 differentiation; as this method considered spatial autocorrelation and is more able to detect low levels of genetic differentiation, it more accurately reflects the true k (33).
Congruent results from our genetic clustering analyses (k = 1) suggest there is no population structuring between individuals in the Central corridor and individuals in the Cockscomb Reserve. The very low levels of genetic differentiation found in this study could be the result of limited dispersal between sampling localities caused by behavioural characteristics of the species, but more likely caused by the presence of significant isolation by distance and a sampling gap between the two regions. Studies conducted with radio telemetry show that jaguars depend on large patches of habitat and can have home ranges that surpass 100 km 2 (3,21,34) however, females have smaller home ranges and tend to avoid roads and human-dominated landscapes at a higher degree, showing preference for intact forests (3, 35). Although having large home ranges and the ability to move considerable distances, jaguars tend to avoid human-dominated areas and show gender-specific differences (3, 35). Genetic subdivision across the country has been

Local-scale connectivity
The relatively high levels of gene flow and low genetic differentiation found in our study attest to the success of the corridor established to connect these two areas of Belize, which were a continuum of jaguar habitat in the distant past. These results are especially informative to aid conservation efforts in other areas of the species range, such as those of the Atlantic Forest of South America, where there is a lack of genetic connectivity among isolated remnant jaguar populations (36). However, anthropogenic barriers (such as Hummingbird highway) could be altering gene flow between core jaguar areas and should alert conservation managers to improve connectivity in future conservation actions and corridor management. The negative impact of roads on jaguar populations should be especially taken into consideration to improve existing corridors or design new ones; road construction and/or expansion within protected areas increases the accessibility of hunters to jaguars and their natural prey, and leads reduction of the potential of using these lands to sustain viable populations of top predators (37). Rabinowitz & Zeller (2010) conducted a range-wide model of landscape connectivity to identify potential corridors that connect jaguar populations across the Americas. Their study provided critical information for conservation actions such as corridor design across the range of the species. However, even if extremely useful for large-scale planning, the model proves more challenging for local or regional corridor design and zoning. The model depends on a least-cost path analysis that relies on coarse-grain environmental data to determine habitat connectivity and could ignore factors that affect how animals utilise the landscape (38). This range-wide corridor could be improved with fine-scale studies that advise targeted-conservation actions (12,38).  Figure 1A). Currently, jaguars seem to move across the two sites, but this highway, other roads and urbanisation, in general, could be shaping population structure by presenting physical barriers to gene flow. To improve connectivity between these sites, the coverage of the Central Corridor needs to expand so that its boundaries cover the areas that would more easily facilitate jaguar movement as evidenced in this study. Conservation efforts should focus on habitat restoration of corridor networks that increase the resistance surface linking Cockscomb Reserve and the Central Corridor to secure movement between and across jaguar core areas could include building wildlife crossings where the resistance surface for movement breaks (e.g. highway junctions in Cayo District and Stann Creek District). The uncertainty over the dispersal ability of jaguars and the extent of use of corridors highlight the importance of incorporating data at a regional scale to better delineate corridors that facilitate gene flow (9,39).

Conclusion
Our results provide a screenshot of genetic patterns of animals whose scats were sampled Linkage disequilibrium (LD) between pairs of loci was performed in GENEPOP v. 4.2 (50) with default settings. We used GenAlex v6.0 (51) to explore multivariate patterns of molecular diversity relative to populations via Principal Coordinate Analysis (PCoA) and Mantel tests of matrix correspondence to test for Isolation by Distance (IBD);; 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 (52), which assigns individuals to a number K of genetically homogeneous groups, based on the Bayesian estimate in accordance to the expected Hardy-Weinberg equilibrium and absence of linkage disequilibrium between loci. We ran STRUCTURE with the LOCPRIOR option 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 (30) using POPHELPER package in R v.1.2.1 (53). Furthermore, we inferred spatial genetic structure with TESS 2.3.1 (54). 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 (55). We tested the CAR, and BYM models with linear trend 16 surface to define the spatial prior for admixture (31); 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 used GENELAND v4.0.6 (55) 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 the genetic diversity and structure among individuals, we reduced the dimensions via a Discriminant Principal Component Analysis (DAPC) without a priori group assignment using the ADEGENT package in R. v2.0.1 (49,56). 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 (57). This multivariate method is useful to identify clusters of genetically related individuals when group priors are lacking. Estimation of clusters is 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 (58). Pairwise relatedness was calculated using the estimators described by Queller and Goodnight, 1989;Li et al., 1993;Ritland, 1996;Lynch andRitland, 1999 andWang, 2002), as well as the dyadic likelihood estimator described in Milligan (2003)

Landscape Connectivity
We predicted the most effective corridor via least-cost path analysis using GDISTANCE package in R v1.1 (66). This approach offers the shortest cost-weighted distance between two sampling points; the program allows calculating grid-based distances and routes and is comparable to ArcGIS Spatial Analyst (67), GRASS GIS (GRASS Development Team 2012), and CIRCUITSCAPE (68). The package implements measures to model dispersal histories and contains specific functionality for geographical genetic analyses (66). The least-cost path analysis was inferred from the most northern GPS point in the Central Corridor to the most southern point in the Cockscomb Reserve. Additionally, we used CIRCUITSCAPE v3.5 (68) to model resistance surfaces of the landscape as an alternative to the least-cost path analysis, which assumes that gene flow is associated to total cost along a single optimal path (69). Circuit theory considers all possible paths and is useful to assess different interactions between different landscape features. With these two approaches, we were able to identify the most probable routes for dispersal and gene flow between localities.
Ecosystem preference costs were based on literature review and expert opinion of habitat use by jaguars (3, 5,6,35,[70][71][72]. To model each ecosystem as analogous to an electrical circuit, each pixel was assigned a resistance value in a scale of 0-9 based upon land cover ( Table 2). The resistance value represents the relative effort required to move from one point to another, and the map of resistance values is used to derive all possible pathways for jaguar movement. Spatial data were obtained from the Biodiversity and  Table 2. Ecosystem types and cost values for jaguar movement based on expert knowledge. Values range from 0 (highly costly) to 9 (no cost for movement).