Fine-scale landscape genetics unveiling contemporary asymmetric movement of red panda (Ailurus fulgens) in Kangchenjunga landscape, India

Wildlife management in rapid changing landscapes requires critical planning through cross cutting networks, and understanding of landscape features, often affected by the anthropogenic activities. The present study demonstrates fine-scale spatial patterns of genetic variation and contemporary gene flow of red panda (Ailurus fulgens) populations with respect to landscape connectivity in Kangchenjunga Landscape (KL), India. The study found about 1,309.54 km2 area suitable for red panda in KL—India, of which 62.21% area fell under the Protected Area network. We identified 24 unique individuals from 234 feces collected at nine microsatellite loci. The spatially explicit and non-explicit Bayesian clustering algorithms evident to exhibit population structuring and supported red panda populations to exist in meta-population frame work. In concurrence to the habitat suitability and landscape connectivity models, gene flow results supported a contemporary asymmetric movement of red panda by connecting KL—India in a crescent arc. We demonstrate the structural-operational connectivity of corridors in KL—India that facilitated red panda movement in the past. We also seek for cooperation in Nepal, Bhutan and China to aid in preparing for a comprehensive monitoring plan for the long-term conservation and management of red panda in trans-boundary landscapes.


Results
Forty-eight candidate maxent models were generated (Table S3) and model with lowest Delta AICc informative AUC (0.911 ± 0.098) was selected, which predicted species distribution better than the random model (Table S1). The generated binary map discriminated between suitable and unsuitable habitat of red panda based on 10-percentile training presence thresholds (Fig. 1b). The jackknife test revealed that bioclimatic variable siwb_bio19, precipitation of coldest quarter contributed 45.6% and maximum temperature of the warmest month contributed 29.4% to predict habitat suitability. The land use land cover (LULC) contributed 12.9% followed by canopy height 11.8% and isothermality found to be less important and contributed only 0.3% (Table 1). The Jackknife test for regularized training gain in the present model showed that the variable with the highest gain, when used in isolation, was precipitation of coldest quarter, which therefore appeared to have the most useful information by itself. The environmental variable that decreased the gain the most when it was omitted, was Max Temperature of Warmest Month, which therefore appeared to have the most information that was not present in the other variables (Fig. S1).
Habitat suitability and landscape connectivity. The study covered 10 protected areas (PAs) that  Table 2). The habitat suitability demonstrated that from West to East, SNP was connected with BRS and BRS was connected with KNP which is the only National Park in Sikkim. The MWLS in south Sikkim was connected with KNP on its northern border. On eastern part of the landscape NVNP and PWLS were connected, however, large amount of habitat suitability was observed outside PA on the eastern part of Sikkim.
The ensemble approach by combining both the genetic divergence and environmental conductance depicted the landscape connectivity and its functionality to support movement corridors of red panda in KL-India. The predicted genetic divergence among the identified clusters was in coherence to the habitat suitability model (Fig. 1b,c). The landscape genetics model by an ensemble conductance surface, lumping both the genetic divergence and environmental suitability demonstrated high-density currents radiated from SNP and NVNP in the North West Bengal (Fig. 1d). Further, to decipher the direction of gene flow among the clusters, we estimated the contemporary migration rate using BayesAss (Fig. S4). The landscape connectivity model and the contemporary migration rate, jointly inferred the permeability of floating currents, indicating structural connectivity in KL-India (Fig. 1d), well supported by the asymmetric gene flow from West to East (Fig. S4). The landscape connectivity model did not show horizontal connection from West to East in KL-India. Instead, the current flow forming a crescent arc to connect KL-India from the West to East with high-density currents radiating from SNP and NVNP of North West Bengal. The current formed bottleneck position of corridor in the KNP of north Sikkim. On the eastern part of the KL, the cumulative current flow was narrow running in longitudinal extent between PAs (Fig. 1d).
Individual identification and assessment of genetic variability. Allelic Drop out (ADO) ranged from 0 to 0.294 while we did not observe any significant amount of False allele (FA) (  www.nature.com/scientificreports/ CRP385 and CRP367 deviated significantly from HWE and only one pair of loci out of 36 pair-wise comparisons were in significant Linkage disequilibrium (LD) (P < 0.0012) after Bonferroni correction. HWE deviation of these markers might be due to the multiple reasons, like inflated frequencies of null alleles, Wahlund effect, and inbreeding due to consanguineous and assortative mating 31 . The mean F IS was 0.276, indicating a significant inbreeding and loss of heterozygosity (Table 3).
Population genetic structure. The STRU CTU RE analysis detected three clusters based on the mean ln P (k) and Delta K (Fig. 2, Fig. S3). We showed probabilities of population assignment at varying K ( Fig. S3). At K-2, a majority of individuals from WS (West Sikkim) showed a distinct cluster but no geographical assignment was observed among the individuals originated from SNP, ES (East Sikkim), and NVNP. At K-3 individuals from SNP showed a distinct isolated cluster while individuals from ES and NVNP assigned to one another cluster. However, a further increase in the K did not reveal any structuring in population. The results indicated that NVNP and ES clusters were relatively connected meta-populations as found shared ancestry in one-another populations while the other two clusters SNP and WS shared ancestry with additional genetic influxes of an unknown population which plausibly be the North Sikkim (NS) and or the historic/uncaptured gene flow from Nepal, which could not be addressed due to unavailability of samples ( Fig. 2a). In congruity, GENELAND also testified to the similar clustering patterns, NVNP and ES grouped into a single cluster and SNP and WS in a distinct clusters (Fig. 2b). The sPCA revealed a west-east differentiation from SNP to NVNP in the allele frequencies, indicating samples of NVNP and ES in one cluster. SNP in the other and WS was assigned as an intermediate population with more affinity to SNP (Fig. 2c, Fig. S2). The DAPC also identified three major clusters, providing strong signals to support red panda populations to exist in metapopulations framework (Fig. 2d). The non-Bayesian methods of population assignment in coherence to Bayesian clustering methods supported strong population genetic structure with asymmetric gene flow among the habitat patches in KL-India.

Gene flow and detection of migrants.
In AMOVA, we found relatively lower genetic variance within population which indicated gene flow within population was higher than the between groups. Red panda populations in KL-India exhibited 58% variance within population than 12% between populations and 30% among group with high Fst (0.321). Though Fst value was not statistically significant but observed signal was an indication of the red panda population to exist in metapopulation framework. Pairwise F ST based gene flow revealed SNP and NVNP populations were highly differentiated (F ST 0.352) and NVNP and ES did not qualify to be separate populations (F ST − 0.052) ( Fig S4). The contemporary migration assessed by the BAYESASS resulted to favour asymmetric movement of red panda from West to East in KL-India (Fig. S4). We detected significant asymmetric migration of individuals from SNP to WS (5.55%) and WS to ES (5.48%) but did not obtain a rational backward gene flow. In contrary, we detected significant bidirectional migration between ES to NVNP (5.91%) and NVNP to ES (6.05%). www.nature.com/scientificreports/

Discussion
Habitat suitability and corridor connectivity. The red panda population in KL has experienced massive habitat loss and fragmentation in the past 11,12,16,32 which must have deteriorated connecting corridors in this landscape. For habitat suitability, predictors such as precipitation, temperature and vegetation with high weightage appeared reasonably correct as vegetation cover, and climatic factors are directly linked to the species diet, survival and reproductive necessities. A few earlier studies addressed the species distribution modelling (SDM) of Red panda in Himalaya 16,20,[33][34][35] . These studies considered different variables like elevation, slope, aspect and distance to water etc. for predicting red panda habitat. These variables may facilitate in old growth forest dominated by Betula utilis, Rhododendron spp. and Abies spp. with dense bamboo cover in the understory and high densities of fallen logs and tree stumps at ground level. Further, the growth of bamboo understory in the temperate forests is highly influenced by the rainfall and temperature 36 and red pandas as being the canopydwelling species prefer temperate forests in the central and eastern Himalayas with dense bamboo undergrowth. Therefore, rainfall, temperature and vegetation cover played a significant role in predicting red panda habitat. Earlier, Choudhury 11 estimated 1,700 km 2 potential habitat of red panda in Sikkim and 200 km 2 in North West Bengal by undertaking exploratory field surveys based on one time efforts. Ziegler 37 estimated 684 km 2 suitable forest for red panda in Sikkim with limited efforts confining to PWLS only. However, Ghose and Dutta 38 estimated 1,200 km 2 suitable forest for red panda in Sikkim and 250 km 2 in North West Bengal. These authors collected primary data by undertaking field surveys through transect and trail monitoring 39 from Sikkim and North West Bengal and also recorded indirect evidences such as droppings, browsing marks, nest sites, pugmarks, skins or pelts, etc. Ghose and Dutta 38 also collected secondary data through questionnaires and public interviews in the villages around the potential red panda habitats. Our field surveys and efforts were relatively comparable to Ghose and Dutta 38 and we obtained an estimated suitable habitat 1,097.26 km 2 in Sikkim and 212.28 km 2 in North West Bengal, which represented a relative decline the habitat suitability in last ten years when compare to Ghose and Dutta 38 . In northern region, red panda habitat was relatively fragmented and circuitscape also supported lack of horizontal connectivity in KL-India. However, the identified corridor favoured red panda movement from West to East in KL India, forming a crescent arc (Fig. S4). The model demonstrated that the identified corridor in the east, south Sikkim and eastern part of north Sikkim was relatively narrow and passing through the non-PAs (Fig. 1d). We propose corridor passing through the non-PAs to be monitored for the emerged developmental activities and enriched by bamboo plantation to avoid further loss of red panda habitat in KL-India. Population genetic structure and gene flow. Explicit Bayesian and non-Bayesian clustering methods to a great extent, showed similar patterns clustering red panda into at least three populations/meta-populations in KL-India. Indisputably, the NVNP and ES was found to represent one single population, However, SNP was another distinct cluster where WS was assigned as an intermediate population with more affinity to SNP. SNP and NVNP that lacks horizontal connection, were two distinct clusters (F ST = 0.352). These populations earlier reported to exist in relatively high density (1 individual per 1.67 sq. km in SNP 20 ; 32 animals in SNP and 34 in NVNP 42 ). Individuals from these populations showed asymmetric gene flow and contributed migration to Sikkim population from either-side. The WS population was rationally distinct from SNP (F ST 0.142), however ES and NVNP, to a large extent represent a single population based on the multiple evidences, e.g. a negative F ST − 0.052 and a high migration rate from one another population (ES to NVNP-5.91% and NVNP to ES-6.05%; Fig. S4). In concurrence to this, AMOVA also supported a high variance (within group) in red panda populations. The observed scenario of contemporary gene flow patterns in KL-India suggested that corridors connecting these populations were functional in the recent past. The pragmatic asymmetric migration detected in the last two generations exhibited asymmetric migration of individuals in KL-India from West to East, SNP to WS and WS to ES. The ES and NVNP had bidirectional migration, also an indicator to support them for being a single population. The observed asymmetric gene flow among populations was due to the landscape heterogeneity and habitat suitability. The circuitscape results evident to prove that SNP and NVNP were two relatively dense populations in KL-India, where individuals might be moving from more stable, relatively high-density population to neighbouring low-density populations 43 . The high migration rate from SNP and NVNP maintained the demographic connectivity in the landscape. Both, the significant genetic differentiation in some areas and the presence of a weak population structure in other regions in KL-India may be explained by the corridor connectivity model (Fig. 1d), which coincides to the results of BAYESASS and supports current flow from West to East forming a crescent arc in the KL-India.
Functional connectivity and conservation priorities. The present study unveils the facts that landscape features have shaped the current distribution patterns of red panda in KL-India. The results demonstrated that NVNP and PWLS were connecting North West Bengal-Sikkim border. But there has been a large gap between PWL to KAS which extended up to the SRS. Further, the suitable habitat was patchy and connecting corridors were narrow between those habitat patches. Habitat on the eastern part of the study landscape was www.nature.com/scientificreports/ more vulnerable to destruction due to the anthropogenic activities. The FLWLS had only a small suitable area with poor connectivity to the other populations. Ziegler 37 conducted a survey in FLWLS but failed to confirm the presence of red panda. There have been recent historic records of red panda from SWLS in West Bengal 44 but intensive surveys need to be conducted to assess the current status of red panda in SWLS, which is also poorly connected to SNP. Though, there are several PAs on the eastern part of KL-India, there is a need to extend PA network in East Sikkim, South Sikkim and eastern part of North Sikkim. Moreover, the territorial forest divisions must take adequate measures, as much of suitable habitat of red panda exists outside the PA. The present study has laid down the foundation to extend the PAs boundaries and by effectively implementing the eco-sensitive zone planning adopted by the Ministry of Environment and Forest, 2011 of Government of India 45 . We propose buffer zone may be declared around the PAs and community conservation area to protect important wildlife corridors. Since, the red panda across the range has experienced habitat loss, fragmentation and population decline due to change in the land use pattern and anthropogenic activities in past few decades. Any significant change in the climatic isotherm might result in vacating the site and or shifting the species to other sites based on varying extent of species resilience and inherent adaptive plasticity 46,47 . Thus, red panda being ecological specialist, serves a good model to test the composite impact of landscapes, historical climate change and contemporary human activities on the possible shift in the ranges.
Further, studying wildlife using genome-wide markers (e.g. GWAS- 15,48 ; SNPs- 15,49 ) is fascinating to evaluate fine scale population genetic structure and investigating loci under natural selection facilitating populations to adapt in the changing climatic conditions 50 . However, to check the immediate effects of landscape features on the genetic variability and population contiguity, the assessment of wild populations using microsatellites is still most cost effective and widely applied way to genetic monitor of free ranging populations 30,51 . Further, consolidating landscape connectivity through mapping corridors, validating movement through genetics and expanding natural PAs is fundamental to make red panda conservation a long-term success across the distribution range. We also seek for the possible collaboration in Nepal, Bhutan and China to aid in preparing a comprehensive monitoring plan for red panda conservation in TBL and to evaluate the ongoing conservation efforts across the political boundaries.

Study design and occurrence data.
Firstly, we stratified the study area based on the forest types, topography to cover all the logistically possible habitat patches in KL-India. Secondly, we adopted a landscape approach, instead of prioritizing PAs for sampling and we aimed to cover maximum reported habitats of red panda in KL-India. Thirdly, for recording the species presence, we adopted a three-pronged approach i.e. transects/trail surveys, questionnaire surveys, and also the remote camera traps following Pradhan et al. 20 , Buckland et al. 39 and Joshi et al. 54 . Further, it was not logistically feasible to cover the entire area systematically due to the rugged terrain, unpredictable weather and poor resources availability. Hence, representative sampling was carried out in KL-India that represented 109 spatially distinct locations, 12 locations from camera traps, seven direct sightings 12 presence locations from reliable interviews and 78 distinct sites where we collected 234 faecal samples of red panda. For model building, we used 56 spatially independent locations after correlation testing among the location following Kramer-Schadt 55 . Informed consent was obtained from all participants who responded to the questionnaire interviews and all methods including human participation in the study was in accordance to the relevant guidelines and regulations of Zoological Survey of India (ZSI), Kolkata. All experimental protocols and methods were approved by the Research and Academic committee (RAC) of ZSI, Kolkata that has a ethics committee/institutional review board. All field surveys and sampling were conducted with taking necessary permits issued by the Forest Department of West Bengal and Sikkim.
Ecological modelling. Selection of environmental variables. Species are sensitive to habitats as well as the climatic isotherms 56 while, any significant change in climatic isotherm to a species having special requirements might result in local extirpation or shifting the species to other ranges 46,47 . Therefore, it is imperative to use variables for SDM which define the likely habitat of a species based on the field observation. Considering these facts, we selected 26 variables out of which 19 variables represented climatic isotherm, four-topography and three variables represented land covers. These variables were also related to the habitat selection of red panda 16,35,57 . The selected variables represented the present environmental conditions and habitat covariates to facilitate red panda distribution patterns in KL-India. The biotic predictors (19) were obtained from the WorldClim data base at 30 arc second scale 58 ; https ://www.world clim.org/ (Table S2)  www.nature.com/scientificreports/ elevation model. For obtaining the forest canopy classes, the satellite data (Landsat 8) was downloaded from US Geological Survey (https ://earth explo rer.usgs.gov/). The forest canopy was classified into four categories viz., very dense, moderately dense, open forest and scrubland 59 . Moreover, considering the fact that the species is a canopy dweller 11 the Canopy height data was obtained from Oak Ridge National Laboratory (https ://webma p.ornl.gov/ogc) for understanding the influence on the species habitat suitability. Since, different variables like climate, landscape, topography and anthropogenic influences used in present study were compiled from various resources present in different resolution scales. Hence, to bring uniformity in the selected data, we re-sampled all raster on the same resolution scale (1 km resolution).
Habitat modelling and landscape connectivity. We employed maximum entropy algorithm for species distribution modelling (MaxEnt 3.4.1), using presence data only [60][61][62] . Bias file was prepared in ENMeval; R package 63 following the methodology suggested by Dudik et al. 64 . We omitted highly correlated variables with a threshold of 0.8 65,66 . Although we selected initially 26 variables but final model was predicted using 19 variables which passed on the multicollinearity test 67 and other related variables with a threshold of 0.8 were excluded to get rid of the over fit model (Table S6) 66,68 . For modelling, we have used 70% of the total locations as training and remaining 30% for testing the significance 60,66,69 . The Akaike's information criteria (AICc) values were used to determine the best fit model with the lowest values in ENMeval. The resultant habitat was classified into suitable, with probabilities 0.584-1 and unsuitable below 0.584 based on the 10-percentile training presence thresholds following Radosavljevic and Anderson 70 (Table S4). The Jackknife was used to analyse the contribution of each variable provided to the MaxEnt model. We tested the model with the receiver operated characteristics area under curve (AUC) where, AUC value > 0.75 indicated high discrimination performance 71 . We also estimated the true skill statistic (TSS-0.742) that compensates for the shortcomings of kappa while keeping all of its advantages following Allouche et al. 72 We generated a raster surface of genetic divergence using Genetic landscape GIS Toolbox Ver. 10.1 73 by generating the inverse distance weighted interpolation within the study area boundary 74 . To evaluate the operational-structural connectivity by predicted movement corridors of red panda in KL-India, we used an ensemble approach by combining both the genetic divergence and environmental conductance (SDM output) in the final circuit model following Mateo-Sánchez et al. 75 and Roffler et al. 76 in Circuitscape 4.0 77 .

Population genetic analyses. DNA extraction, PCR amplification and microsatellite genotyping.
We collected 234 red panda feces, i.e. 87 feces from North West Bengal and 147 from Sikkim in KL-India (Fig. 1a). All samples were stored in 70% ethanol and DNA was extracted using QIAamp Fast DNA Stool Mini Kit (QIAGEN Germany) following manufacturer's instructions. Nine polymorphic STRs, i.e. Aifu01, Aifu05 from Liang 78 and CRP357, CRP385, CRP381, CRP367, CRP409, CRP240 and CRP260 from Yang 79 were amplified into two multiplex PCRs. Forward primer of all nine microsatellite loci was fluorescently labelled at 5′ with one of the four dyes, FAM, VIC, NED and PET (Table S5). The PCRs were carried-out in 10 µl reaction volume following QIAGEN Multiplex PCR Kit (Qiagen, Germany). The thermal cycle profile was: initial denaturation at 95 °C for 15 min, followed by 40 cycles of PCR and a final step of 72 °C for 30 min. The annealing temperature (Ta) for multiplex 1 (CRP381, CRP367, CRP409, CRP240 and CRP260) was 55 °C and multiplex 2 (Aifu01, Aifu05, CRP357, CRP385) was 57 °C (Table S5). The PCR products were resolved on an ABI 3730 Genetic analyzer (Applied Biosystems, Foster City, CA, USA) and allele scoring was done using Gene Mapper 4.1 (Applied Biosystems, Foster City, CA, USA).
Genotyping error and individual identification. We genotyped each sample four times to minimize genotyping errors and a heterozygote was ascertained only if there were different alleles in at least three independent attempts. In addition, we also followed an independent allele scoring method where two different researchers perform allele scoring individually and only the consensus genotypes were used for further analysis. The genotyping errors arising due to null allele and the presence of stutters, scoring errors were assessed using MICRO CHECKER 2.2.2 80 . Maximum likelihood allele dropout (ADO) and false allele (FA) error rates were quantified using PEDANT version 1.0 involving 10,000 search steps for enumeration of per allele error rates 81 . To avoid ambiguity in ascertaining unique genotypes, we limited the number of loci used based on their high success rate (> 90%), presence of no or minimum genotyping errors and exhibiting an informative P ID value (probability of obtaining identical genotypes between two samples by chance). The locus wise and cumulative probability of identity for unrelated individuals (P ID ) and siblings (P ID sibs) were calculated following identity analysis module in GenAlEx version 6 82 .
Genetic diversity and inbreeding. The genetic diversity estimates were accounted by calculating the number of allele (Na), observed (Ho) and expected (He) heterozygosity using GENALEX 6 82 . For the Hardy-Weinberg equilibrium (HWE) test, we followed the probability test approach using the program GENEPOP version 4.2.1 83 . Wright's inbreeding coefficient (F IS ) was estimated following Weir & Cockerham 84 using GENEPOP 83 . Linkage disequilibrium (LD) was tested using GENEPOP 83 to determine the extent of distortion from independent segregation of loci following 10,000 dememorizations, 500 batches and 10,000 iterations per batch after Bonferroni correction 85 . www.nature.com/scientificreports/ different Bayesian clustering methods, individuals were assigned exclusively on the basis of their multi-locus genotypes (e.g. STRU CTU RE), and also using the both, multi-locus genotypes and geo-referenced information (e.g. GENELAND). Non-Bayesian multivariate ordination analyses, i.e. discriminant analysis of principle components (DAPC) and spatial principle component analysis (sPCA) were also used, to compare population assignment with the Bayesian clustering outputs 87 since they were not based on any model assumptions.
In Bayesian analysis, STRU CTU RE 2.3.4 software 88 was used to determine the number of genetic clusters (K) following 20 iterations (20,000 burn-in; 200,000 Markov chain Monte Carlo replicates in each run) with NOPRIOR with admixed and correlated allele frequencies. We considered there were K populations (1 to 10), with repeating each analysis for 10 times at each K value. The most probable cluster was calculated via estimating the distribution of Delta K 89 using STRU CTU RE HARVESTER v.0.68 90 . GENELAND v 4.0.3 91 was run through an extension of R v.3.0.1 with the correlated allele frequency and spatial uncertainty model. We allowed K to vary between 1 and 10 following 20 independent runs, each with 100,000 iterations, and a thinning of 1,000. The inferred spatial clusters were georeferenced in ArcGIS 10.6. Two non-Bayesian clustering methods i.e., DAPC and sPCA, were run to assign the possible clusters in Adegenet v1.3.4 package of R 92 .
Gene flow and migration rate. We assessed Fst based gene flow in last 150-200 years 93 among the different sub-populations of red panda in KL-India using Arlequin v 3.5.2 86 . The contemporary and asymmetric migration rate in last two generations were estimated using BayesAss 1.3 43 . We used 9 × 10 6 iterations, with a burn-in of 10 6 iterations, 1,000 number of permutations and a sampling frequency of 2000 to ensure that the model's starting parameters were sufficiently randomized. We also estimated the first-generation migrants between all pairs of subpopulations of red panda in KL-India using GENECLASS 2.0 94 , an approach described by Paetkau 95 for likelihood computation (Lhome/Lmax), with 1,000 simulations at an assignment threshold (alpha) of 0.01 and 0.05.

Data availability
The multilocus genotype data of 24 unique individuals is available on DRYAD (https ://doi.org/10.5061/dryad .2280g b5nz) and all other relevant data is included in the manuscript. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.