Anthropogenic pressures drive population genetic structuring across a Critically Endangered lemur species range

In recent decades Madagascar has experienced significant habitat loss and modification, with minimal understanding of how human land use practices have impacted the evolution of its flora and fauna. In light of ongoing and intensifying anthropogenic pressures, we seek new insight into mechanisms driving genetic variability on this island, using a Critically Endangered lemur species, the black-and-white ruffed lemur (Varecia variegata), as a test case. Here, we examine the relative influence of natural and anthropogenic landscape features that we predict will impose barriers to dispersal and promote genetic structuring across the species range. Using circuit theory, we model functional connectivity among 18 sampling localities using population-based genetic distance (FST). We optimized resistance surfaces using genetic algorithms and assessed their performance using maximum-likelihood population-effects mixed models. The best supported resistance model was a composite surface that included two anthropogenic features, habitat cover and distance to villages, suggesting that rapid land cover modification by humans has driven change in the genetic structure of wild lemurs. Primary conservation priority should be placed on mitigating further forest loss and connecting regions identified as having low dispersal potential to prevent further loss of genetic diversity and promote the survival of other moist forest specialists.

Remote sensing and landscape feature selection. We selected five landscape features hypothesized to influence ruffed lemur movement and gene flow: (1) rivers, (2) topography, (3) roads, (4) habitat type, and (5) proximity to human habitation (Fig. 2). Categorical surfaces analyzed included habitat type (primary forest, degraded habitat, and matrix), rivers, and roads, while continuous surfaces included topographic position index (TPI), and distance to nearest village (in meters). Habitat type was derived by reclassifying a vegetation raster from The CEPF Madagascar Vegetation Mapping Project (http://www.vegmad.org/) to three primary habitat types such that primary forest represented a combination of humid and littoral forests, degraded habitat was derived from the degraded forest classification, and matrix was derived from all other categories. River and road data were downloaded from DIVA-GIS (http://www.diva-gis.org/), clipped to include only major features (i.e., using only roads that had been designated as 'roads' rather than 'trails'), and then rasterized. Topographic position index, a measure that compares the elevation of each cell in the raster to the adjacent landscape and calculates a quantitative value that is indicative of the cell's relative position (i.e., slope, valley, plain, or ridge), was generated from a 30 arcsecond resolution digital elevation model (DEM; downloaded from DIVA-GIS at http://www. diva-gis.org/) using the Topography Tools toolbox in ArcGIS v10.3.1 78 . Village locations were obtained by selecting point locations with a designation of "Populated Places" from a gazetteer of foreign geographic feature names (data was downloaded from DIVA-GIS at http://www.diva-gis.org/gdata; source: GEOnet Names Server at the U.S. National Geospatial-Intelligence Agency). Village proximity was derived by calculating the distance of each cell in a 1 km resolution raster to the nearest village center. All surfaces were converted to a uniform geographic coordinate system (Universal Transverse Mercator; UTM), resampled at a resolution of 1 km, and clipped to the study extent for analysis. Previous work has shown that changes in spatial resolution do not significantly alter the  (left). Colored nodes correspond to current subspecies status, as indicated to the right of the map. We provide results from an earlier structure plot (right) modified from Baden et al. (2014), which identify the Mangoro River as a likely driver of V. variegata population genetic structure into northern (red) and southern (green) genetic clusters.
www.nature.com/scientificreports www.nature.com/scientificreports/ results of landscape genetic analyses 74,79 . Therefore, this resolution was chosen as a trade-off between retaining detail across the landscape and minimizing processing time for analyses.
Distance is implicitly incorporated into the resistance distances calculated by Circuitscape, and thus straightline Euclidean distance between sites (i.e., geographic distance) was not included as an additional factor in our models. We did, however, conduct a linear regression between pairwise Euclidean and genetic (F ST ) distances to identify the proportion of our data influenced by geographic distance alone.
Landscape resistance parameterization. Resistances between sampling localities were calculated in Circuitscape 4.0 77 based on average pairwise resistances using an eight-neighbor connectivity scheme and optimized using the R package ResistanceGA 39 . ResistanceGA utilizes genetic algorithms to adaptively search a broad parameter space to determine the optimal resistance values that best describe pairwise genetic differentiation (in our study, F ST ). Inspired by principles of biological evolution, genetic algorithms create a population of individuals with traits (parameters to be optimized) encoded on "chromosomes". The fittest individuals from each generation (i.e., those with genotypes or parameter combinations that solve the fitness function) survive to reproduce 80 . Parameter space is explored via 'mutations, ' whereby new parameter values are generated, as well as via 'crossover' , whereby genetic information is exchanged. The population continues to evolve -that is, parameterization of the resistance surfaces run -until a sufficient number of generations have passed without an improvement in fitness 81 (see ref. 39 for details). This approach makes no a priori assumptions about the direction or magnitude of the resistance between landscape and genetic distances, allowing for a more thorough investigation of the relationship between landscape features and gene flow than more widely used methods (i.e., expert based value assignments 34,82 ).
Continuous surfaces were optimized using Monomolecular and Ricker transformations, while categorical surfaces were optimized by holding one feature constant at a value of 1 and then adjusting resistance values for , while river and road data were downloaded from DIVA-GIS (http://www.diva-gis.org/), clipped to include only major features (i.e., using only roads that had been designated as 'roads' rather than 'trails'), and then rasterized. Distance to nearest village was generated by selecting point locations with a designation of "Populated Places" from a gazetteer of foreign geographic feature names (data was downloaded from DIVA-GIS at http://www. diva-gis.org/gdata; source: GEOnet Names Server at the U.S. National Geospatial-Intelligence Agency), and calculating straightline distances from each site to the nearest village locale. Topographic position index was generated from a 30 arcsecond resolution digital elevation model (DEM; downloaded from DIVA-GIS at http:// www.diva-gis.org/) using the Topography Tools toolkits in ArcGIS v10.

Resistance optimization and model selection.
We evaluated the resistance optimization process for each surface (i.e., landscape feature) using AICc (Akaike's Information Criterion corrected for small/finite population size 83 ), which was determined from linear mixed effects models with MLPE parameterization 84 and evaluated by maximum likelihood in lme4 85 . To account for our uneven sampling design, we conducted bootstrap resampling using 75% of the data (13 sampling locations) to control for bias following Ruiz-Lopez et al. 86 ; these data were randomly selected without replacement and then optimized surfaces were fit to the selected data. Following 10,000 iterations, average rank and average model weight (ω) were determined for each resistance surface, along with the frequency with which a surface was ranked as the top model (π ) in order to address uncertainty in the top model; Burnham and Anderson 87 identify π as the bootstrap equivalent of ω . Following the identification of the top surfaces in isolation, we ran a Spearman's rank correlation between the surfaces to assess degree of correlation. We created and optimized composite surfaces by combining the top models; surfaces that had both a greater selection frequency (π ) than distance alone and were selected more than one percent of the time (π > 0.01) were used to generate composite surfaces. All single and composite surface optimization processes were conducted at least twice as per recommendations by Peterman 39 to ensure convergence on the top model(s). Following optimization, we again conducted a bootstrap model selection using 10,000 iterations and average rank, ω, and π were calculated in order to assess composite models in relation to their component surfaces. Finally, current flow across the landscape was visualized in Circuitscape using the best supported resistance surface.
Landscape surfaces generated during the current study are available in the Zenodo repository, as is all code used to parameterize, optimize, and asses resistance layers [http://doi.org/10.5281/zenodo.3519173].

Results
Landscape genetic analysis. We found a significant signature of isolation-by-Euclidean distance, which explained 28.2% (Pearson's r = 0.531) of the observed population genetic structure (Fig. S1). Despite evidence of IBD, two surfaces -habitat cover and distance to nearest village -were much more strongly associated with genetic differentiation than geographic distance, with the former two variables having average model weights more than 28 times higher than the latter (Tables 1 and S1). Resistance to forest cover was lowest in primary forest habitat, slightly higher in matrix, and highest in degraded landscape (Fig. S2); resistance to human habitation decreased with increasing distance from the nearest village, but increased again after approximately 8 km (Fig. S2). The remaining three surfaces (TPI, roads, and rivers) explained slightly more variation than distance, but were seldom chosen as the top model (less than 0.21% of the time) and thus were not considered in the composite analysis.
To assess the combined effects of habitat cover and distance to the nearest village, we created a composite surface that combined both habitat cover and distance to the nearest village and assessed this against the isolated surfaces, as well as straight-line Euclidean distance alone. As with its component surfaces, resistance in the composite surface decreased as distance from the nearest village increased but then increased again over large distances. Resistance in matrix was lower than, albeit almost equal to, primary forest, with degraded habitat again showing the highest resistance. The composite surface had the greatest average rank, highest average weight (ω), and was chosen as the top model 70% of the time, lending support for the composite model as the best predictor of the genetic data ( Table 2). Our visualization of the composite model (Fig. 3) shows a high degree of current flow through remaining sections of V. variegata's habitat across its range − particularly in the central region, where gene flow is likely the greatest -while also illuminating areas of near complete isolation in the south. Regions that are covered in primary forest and are distant from villages tend to be where current density is the greatest, such that the primary flow of current is located in the contiguous forest corridor stretching through the center of the existing V. variegata range.

Discussion
Recent studies have identified a widespread, latitudinally-structured phylogeographic pattern in several of Madagascar's fauna, including its humid-forest mammals 88 . Several factors may be driving this pattern, including  Table 2. Results from bootstrap selection on optimized linear-mixed effects models on composite surface and its component single surfaces K = number of parameters following continuous surface transformation or number of categories in categorical surfaces; Avg. rank = average model rank following 10,000 bootstrap iterations; ω = average model weight averaged over 10,000 bootstrap iterations, representing the probability that the model is the best of the set; π = proportion of bootstrap iterations in which model was chosen as the top model. www.nature.com/scientificreports www.nature.com/scientificreports/ past episodes of forest contraction and expansion and biogeographic barriers (i.e., rivers or valleys). These landscape features are thought to operate on much deeper time scales than anthropogenic change in terms of their influence on an organism's dispersal abilities, population connectivity, and genetic structure 89 . In an earlier study, we identified genetic clustering on either side of the Mangoro-the largest river in eastern Madagascar-lending support to these findings, and leading us to hypothesize that Madagascar's rivers structure V. variegata genetic diversity 61 . We were therefore surprised to find little support for this hypothesis. Rather, our results suggest that anthropogenic factors are the primary drivers of ruffed lemur genetic structure.
There is some debate surrounding the extent to which rivers present complete barriers to vertebrate movement (e.g., [90][91][92][93][94][95], particularly in montane species where river headwaters are still relatively small 96 . However, rivers have been strongly implicated in the biogeography of Malagasy fauna, particularly primates [66][67][68][69] , and it is unlikely that Varecia are exempt from these patterns. Thus, we propose two alternative explanations for this finding. On the one hand, our results may be a product of methodological constraints, including our inability to filter our waterways landscape layer to include only the largest rivers (such as the Mangoro), as freely available data for this feature do not contain sufficient metadata to differentiate between major and minor waterways. Alternatively, we argue that natural signatures of migration and drift are being swamped by patterns caused by anthropogenic activities.
Madagascar is experiencing exceptional threats to its biodiversity 97 . In particular, deforestation has reduced forest cover by 44% between 1953 and 2014 45 , and is continuing at a pace that is projected to eliminate all of Madagascar's eastern rainforest by 2077 70 . Of the lemurs, ruffed lemurs are particularly sensitive to anthropogenic pressure, and are among the first taxa to disappear in light of habitat disturbance 49,60 . Here, we offer further evidence of this, by presenting results indicating that anthropogenic factors alone are driving current patterns of ruffed lemur genetic structure, a finding which has important consequences for their long-term viability.
Animals are more likely to disperse through landcover that closely reflects the habitat in which they evolved 98 . Our results further support this notion, with ruffed lemur dispersal capacity being at its highest in intact primary forest, followed by matrix and degraded forest. That matrix (e.g., agriculture, grassland) proved to be a more suitable conduit to gene flow was surprising, but this finding is likely driven by the orientation of matrix to our ruffed lemur sampling localities. In our study, very few sampling locations were directly separated by matrix. Rather, with the exception of a few sites (e.g., Kianjavato, Manombo), samples were collected from protected, primary rainforest sites that were separated from the surrounding matrix by intervening degraded forest habitat. For this reason, black-and-white ruffed lemurs (or at least those included in this study) rarely encounter matrix, and it therefore did not impose a significant resistance to their movement and dispersal.
Still, it was surprising that matrix proved more permeable than degraded habitat. It is possible that this is an artefect of the classification system used, and that degraded habitat in this classification scheme represents what might be perceived as 'matrix' to ruffed lemurs (i.e., low quality, small-crowned fruiting trees and shrubs), effectively serving as an additional barrier to movement between patches. There is, however, indirect evidence that individuals may cross highly degraded regions at short distances, such as those that segregate rainforest fragments 52 .
The exact dispersal capacity of black-and-white ruffed lemurs is at present unknown; however, maximum mammalian dispersal distances have been shown to scale isometrically with home range size 99 . The average individual black-and-white ruffed lemur home range is ~15-20 ha 54,55 , giving an inferred migration capacity of no more than 200 km. This indicates a large dispersal capacity that would likely reduce any time lag exhibited by historic landscapes. Furthermore, all nearest neighbor sampling localities were located within this maximum dispersal distance, suggesting a high likelihood that ruffed lemurs reach at least one or more sampling locations during dispersal events. Given the orientation of the existing range for this species, however, the distances between non-neighboring sampling localities (up to 861.74 km) was often beyond the assumed functional dispersal distances of the taxon -a sampling strategy that may have decreased the signature of landscape effects on gene flow during population-level analyses 100 . At present, there are clear limitations in researcher ability to directly observe or test ruffed lemurs' dispersal capacity. Indeed, our analysis is among the first to quantitatively relate landscape features to a measure of dispersal using genetic data. With rapid improvements in telemetry, it is now feasible to track individuals via GPS collars long-term, allowing for the possibility of capturing dispersal events and gathering spatially-explicit movement data moving forward (EEL personal observation).
As with habitat degradation, we also found increased resistance to movement with increasing proximity to human habitation. Curiously, gene flow increased again after approximately 8 km of isolation; this is likely because most villages in Madagascar are situated well within 8 km of forest, such that dispersing ruffed lemurs would rarely encounter habitat located at distances greater than 8 km from human habitation.
Unfortunately, we were unable to test the impacts of historical forest cover on the genetic structure of ruffed lemurs, as landscape maps such as these are not georeferenced and publicly available. This is an important consideration, given that historical habitat loss may yet impact future primate populations due to time lags in population responses 101 . Several factors influence the degree to which observed genetic structure reflects contemporary versus historical landscapes, including generation time, dispersal distance, and population size and demographics of the study species, as well as genetic metrics chosen for analysis 102 . Black-and-white ruffed lemurs have relatively fast generation times for primates of their body size (approximately 7-8 years 52 ), suggesting that an influence from historic landscape features may no longer be detectable in contemporary genetic structure 102 . Moreover, although forest loss was already well-advanced as early as the 1600s, most significant habitat modification has occurred in the past 70 years 44,45 , suggesting that signatures of habitat modification have predominated as drivers of genetic structure in ruffed lemurs in less than ten generations. With vagile species, this generation time has shown to be sufficient for detecting newly introduced barriers to gene flow 103  www.nature.com/scientificreports www.nature.com/scientificreports/ conservation implications. To date, the vast majority of landscape genetic analyses have been conducted in temperate locations 31 , with a significant underrepresentation of tropical regions despite the vast majority of worldwide biodiversity being found therein 41 . Madagascar is considered both a biodiversity hotspot and a conservation priority due to its high levels of endemism in most taxonomic categories for a relatively small geographic area and extremely high loss of primary habitat 41 . Presently, over 90% of endemic lemur taxa are classified as Vulnerable and anthropogenic activities have long been implicated as a primary cause for taxonomic decline 104 . Given their sensitivity to habitat degradation, ruffed lemurs are often used as an indicator species and present an opportunity to assess the influence of anthropogenic activities across a large region of Madagascar. Our investigation suggests that rapid land cover modification by humans (e.g., refs 45,105 ) has driven change in genetic structure of this indicator species, masking any signatures of natural influences.
At present, Madagascar's natural forests cover approximately 8.9 Mha (15% of the national territory) and include 4.4 Mha (50%) of moist forests; however, estimated annual deforestation rates have progressively increased since 2005, and approximately half of Madagascar's forest (46%) is now located at less than 100 m from the forest edge 45 . Connecting regions of low dispersal, particularly those in the southern and northern-most sampling sites, should therefore be targeted as a primary conservation priority in order to connect currently isolated ruffed lemur populations. Prior studies indicate that increased fragmentation and/or small fragment size are related to increases in genetic differentiation and lowered genetic diversity (e.g., refs 52,106,107 ). Higher levels of genetic variability have been described in the northern ruffed lemur populations; however, forest loss and fragmentation have continued in this region and across the country at rates of up to 2.5% habitat lost per year since the collection of these data 105 . Additionally, northern populations experience greater intensities of hunting pressures compared to southern populations 108 , which may drive extreme losses in population size and further reduce genetic variability and increase population differentiation through greater levels of inbreeding or local extinctions. Hunting in southern localities has been recently recognized as a threat to populations as well 109 . Further investigations into the local drivers of genetic structure in this taxon are currently underway and will clarify the regional impacts of natural features and landcover modification on more immediate dispersal capacity, as well as more explicitly assess the impacts of historical landcover of observed genetic structure in this taxon.
Our study uses ruffed lemurs as a test case; however, they also represent a dispersal-, resource-, and area-limited 'umbrella' species (sensu refs 110,111 ) whose habitat requirements for persistence are believed to encapsulate those of an array of associated species. As such, it is likely that the trends we see here may be echoed by other moist forest specialists, either currently or in the near future. One somewhat controversial approach might therefore be to promote ruffed lemurs as conservation 'flagships' (sensu ref. 112 ) for Madagascar's eastern rainforest corridor. Although this concept has faced some criticism (e.g., refs [113][114][115], there is nevertheless evidence of greater 'willingness-to-pay' among private sponsorship programs for conservation focusing on charismatic species 116 . Ruffed lemurs, the largest living members of the lemurid family, can be easily recognized by their "especially long, luxuriant coat relative to other lemurs" 117 and their suspensory feeding acrobatics 118 ; they may therefore be appropriate 'Cinderella (flagship) species' sensu Smith et al. 119 that can be used to attract public support for conservation efforts. Planning conservation action around this taxon, both logistically and for promotional purposes, may be one strategy toward protecting the diverse eastern rainforest species in Madagascar. These and other novel solutions will become increasingly important, as anthropogenic activities including those demonstrated herein continue to overwhelm the impact of natural processes on ecology with far-reaching, long-term and cascading consequences for the whole Earth system 120,121 (reviewed in ref. 122 ).