Terrain Ruggedness and Canopy Height Predict Short-Range Dispersal in the Critically Endangered Black-and-White Ruffed Lemur

Dispersal is a fundamental aspect of primates’ lives and influences both population and community structuring, as well as species evolution. Primates disperse within an environmental context, where both local and intervening environmental factors affect all phases of dispersal. To date, research has primarily focused on how the intervening landscape influences primate dispersal, with few assessing the effects of local habitat characteristics. Here, we use a landscape genetics approach to examine between- and within-site environmental drivers of short-range black-and-white ruffed lemur (Varecia variegata) dispersal in the Ranomafana region of southeastern Madagascar. We identified the most influential drivers of short-range ruffed lemur dispersal as being between-site terrain ruggedness and canopy height, more so than any within-site habitat characteristic evaluated. Our results suggest that ruffed lemurs disperse through the least rugged terrain that enables them to remain within their preferred tall-canopied forest habitat. Furthermore, we noted a scale-dependent environmental effect when comparing our results to earlier landscape characteristics identified as driving long-range ruffed lemur dispersal. We found that forest structure drives short-range dispersal events, whereas forest presence facilitates long-range dispersal and multigenerational gene flow. Together, our findings highlight the importance of retaining high-quality forests and forest continuity to facilitate dispersal and maintain functional connectivity in ruffed lemurs.


Introduction
Animal movement, or an individual's change in spatial location through time, can occur at multiple spatial and temporal scales [1,2]. It is a fundamental aspect of a primate's life and is key to individual survival and fitness, population and community structuring, and species evolution [1,3]. Behaviors involving close-range movement, such as foraging, locating a mate, or avoiding predators, occur periodically at relatively small spatial and temporal scales [1]. By contrast, long-range movement, such as natal or secondary dispersal, typically occurs only once or a few times throughout a primate's lifetime and at relatively cover, deterred by proximity to human settlements, and, contrary to expectations, environmental features such as rivers and altitude are unrelated to range-wide gene flow in the species [46]. It is, however, unclear whether these patterns might vary at different spatial scales (as in [48][49][50][51]). Further work has identified local habitat quality as a major predictor of ruffed lemur occupancy across the species' range [55]; however, studies have not yet assessed the impact of source habitat quality on ruffed lemur dispersal decisions.
Despite evidence of genetic isolation in ruffed lemurs south of the Mangoro River [46], we recently found evidence of functional connectivity throughout Ranomafana National Park and the adjacent Ambositra-Vondrozo Forest Corridor (COFAV), which encompass the greatest contiguous stretch of forested habitat remaining within ruffed lemurs' current southern range [59]. As with all remaining ruffed lemur habitats, these areas are subject to ongoing habitat transformation due to slash-and-burn agriculture (tavy), mining, and selective logging [60][61][62], leading to significant environmental heterogeneity across the region. Furthermore, within ruffed lemurs' habitats, there is evidence of significant environmental variability in terms of forest structure and floristic diversity resulting from historic and contemporary anthropogenic activities [59,62]. This environmental variation across the landscape, combined with documented patterns of ruffed lemur connectivity throughout the southeastern rainforest corridor, makes the Ranomafana National Park and the adjacent COFAV region an excellent location to evaluate the role of environmental variation on patterns of short-range dispersal in this species.
Given ruffed lemurs' reliance on large-canopied trees and high levels of frugivory [63], we expected that higher-quality sites (e.g., those with greater productivity and/or greater structural complexity) would be more attractive to dispersers than sites of lower quality (as identified in [59]). Between-site forest cover is a major driver of dispersal and gene flow in black-and-white ruffed lemurs across the species' range [46]. We, therefore, expected that reductions in forest cover (via loss or degradation) would impede gene flow throughout the Ranomafana region. We tested for the influence of both historic and contemporary forest cover, as time lags are often present when detecting the impacts of environmental variables on genetic signatures [64]. Furthermore, as elevation closely relates to forest structure and floristic diversity worldwide [65][66][67][68][69][70][71], we also expected altitude to indirectly drive gene flow via its effect on plant communities throughout the Ranomafana region. Because environmental impacts on species are often scale-dependent [48][49][50][51], we compared our results to those from Baden et al. [46] to evaluate similarities and differences in environmental drivers of short-range dispersal (this study) and long-range dispersal/multigenerational gene flow in ruffed lemurs, as previously identified [46]. Furthermore, to our knowledge, our study is the first to evaluate the impact of within-site habitat characteristics on gene flow in primates and will inform our understanding of how local environmental characteristics influence the dispersal process, particularly the emigration and immigration phases [4]. Finally, by comparing regional and species-wide drivers of gene flow, this study strengthens our understanding of the role environmental variation plays in community structuring, and by extension, the evolutionary process at varying scales.

Ethics Statement
This research adhered to the American Society of Primatologists Principles for the Ethical Treatment of Non-Human Primates. The research complied with the laws and guidelines set forth by ANGAP/Madagascar National Parks and Hunter College IACUC (#AB-impact 4/18-01).

Landcover Classification and Landscape Feature Selection
We selected eight landscape features that were hypothesized to influence ruffed lemur movement and gene flow: (1) 1990 forest cover (historic); (2) 2016 forest cover (contemporary); (3) rivers; (4) altitude; (5) topographic position index (TPI); (6) terrain ruggedness index (TRI); (7) fire density; and (8) canopy height ( Figure 2). Categorical surfaces analyzed included 1990 and 2016 forest cover and rivers, while continuous surfaces included altitude, TPI, TRI, fire density, and canopy height. Forest cover in 1990 and 2016 were classified using spectral mixture analysis (SMA) and linear spectral unmixing (see Supplementary Methods S1 for details). River data were downloaded from Open Street Map (https://www.openstreetmap.org/; accessed on 5 January 2021) as vector data and were rasterized. Altitude data from the Shuttle Radar Topography Mission (SRTM) were downloaded from the USGS Earth Explorer data portal (https://earthexplorer.usgs.gov/; accessed on 5 January 2021) at 1 arc second (approximately 30 m) resolution. Topographic position index, a measure that compares the elevation of each pixel in the raster to the adjacent landscape and calculates a quantitative value that is indicative of the pixel's relative position (i.e., slope, valley, plain, or ridge), was derived from altitude data using the topographic position index function from the GDAL library (http://www.gdal.org/; accessed on 5 January 2021) in QGIS v.3.20. Terrain ruggedness index, a measure of the ruggedness of a pixel calculated by comparing elevation differences between a pixel and its eight neighboring cells within the raster [81], was also derived from the altitude data using the terrain ruggedness index function from the GDAL library (http://www.gdal.org/; accessed on 5 January 2021) in QGIS v.3.20. Fire data were downloaded from the NASA Fire Information for Resource Management System (FIRMS) from both the MODIS Collection 6 (MC6; 2001-2016) and the VIIRS S-NPP 375 m (2012-2016) datasets. Data from MC6 and VIIRS were subsampled to only include nominal and high-confidence fire reports (above 30% confidence; [82]), data were merged, and fire density was estimated using the Heatmap Kernel Density Estimation function in QGIS v.3.20 with a radius of 500 m (based on suggestions from Page et al. [83]). 'No-data' pixels were filled with a value of 0 to indicate no fires in that pixel. We downloaded 2019 (contemporary) canopy height data from Potapov et al. [84], which were generated through the integration of the NASA Global Ecosystem Dynamics Investigation (GEDI) spaceborne lidar system on the International Space Station and a timeseries of Landsat imagery. We resampled canopy height data to a 10,000 m 2 resolution and corrected non-forest values (over 60) to 0. All surfaces were converted to a uniform geographic coordinate system (Universal Transverse Mercator; UTM), resampled at a resolution of 10,000 m 2 , and clipped to the study extent for analysis. Previous work has shown that changes in spatial resolution do not significantly alter the results of landscape genetic analyses [85,86]. Therefore, this resolution was chosen as a trade-off between retaining detail across the landscape and minimizing processing time for analyses. Straight-line Euclidean distance between individuals (i.e., geographic distance) was not included as an additional factor in our models, as distance is incorporated as a null model in ResistanceGA. We did, however, conduct a linear regression between pairwise Euclidean and genetic (Ar) distances to identify the proportion of our data influenced by geographic distance alone.

Resistance Surface Parameterization and Optimization
Resistances between sampling localities were calculated using the commuteDistance function from the R package gdistance [87], based on average pairwise resistances using an eight-neighbor connectivity scheme, and optimized using the R package Resis-tanceGA [88]. 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, Ar). 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; [89,90]). Con-tinuous 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 all other features between the values of 0 and 5000 for single surfaces and 10,000 for composite surfaces. The analyses were run in parallel on the Extreme Science and Engineering Discovery Environment (XSEDE) Bridges-2 Regular Memory platform [91] using a modified Singularity environment from Finnish Center for Science (CSCfi; https://github.com/CSCfi/singularity-recipes; accessed on 11 March 2021). We evaluated the resistance optimization process for each surface (i.e., landscape feature) using log-likelihood (corroborated with AICc; Akaike's Information Criterion corrected for small/finite population size; [92]), which was determined from linear mixed-effects models with MLPE parameterization [93] and evaluated by maximum likelihood in lme4 [94].

Resistance Surface Parameterization and Optimization
Resistances between sampling localities were calculated using the commuteDistance

Resistance Surface Model Selection
To account for our uneven sampling design and control for bias, we conducted bootstrap resampling using 75% of the data (47 sampling locations) following Ruiz-Lopez et al. [47]. These data were randomly selected without replacement, and then optimized surfaces were fit to the selected data. Following 10,000 iterations, the 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 & Anderson [95] identifyπ as the bootstrap equivalent of ω. After identifying the top surfaces in isolation, we ran Spearman's rank correlation between the surfaces to assess the 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 (π > 1.00) were used to generate composite surfaces. All single and composite surface optimization processes were conducted at least twice, as per recommendations by Peterman [96], 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 to assess composite models in relation to their component surfaces. Finally, current flow across the landscape was visualized in Circuitscape v4.0.5 using the best supported resistance surface(s). Landscape surfaces generated during the current study are available on GitHub (https://github.com/amandamancini/ruffed_lemur_landscape_ genetics accessed on 5 January 2021), along with the Singularity image and all code used to parameterize, optimize, and assess resistance layers.

Gravity Model
To evaluate how within-site conditions influence dispersal and functional population connectivity, we used singularly constrained gravity models based on a saturated network. Gravity models use a network-based approach composed of nodes and edges to evaluate both within-and between-site environmental drivers of functional population connectivity, respectively [97]. To evaluate within-and between-site drivers, gravity models integrate three parameters: the distance between sites (ω), influence of within-site conditions on attraction of individuals to or from a site (υ), and the resistance of intervening landscape features between sites (c). We used proportion of shared alleles (D PS ) as our measure of gene flow as this metric is free of equilibrium assumptions and can represent shared information, therefore satisfying the "maximum entropy information minimizing" approach used in landscape genetic gravity models (for more details see Murphy et al., [98]). We calculated D PS using the R package adegenet and implemented the singularly constrained gravity models in R in package GeNetIt [37].
We used Euclidean distance between dyads (ω) and landscape resistance (c) calculated from optimized resistance assignments in ResistanceGA [96]. Within-site (υ) variables included forest structure and floristic diversity measures, average Normalized Difference Vegetation Index (NDVI), and topographic variation. Forest structure measures included average canopy height, basal area, and stem density at each sampling site, and floristic diversity included effective number of species (ENS) and average Importance Value Index (IVI) of ruffed lemur food trees at each site. Both forest structure and floristic diversity were quantified from three to five 25 m by 25 m plots at each site (described in Mancini [59]). Normalized Difference Vegetation Index (NDVI) was used to evaluate site vegetation density and was calculated as (NIR − VIS)/(NIR + VIS) from the 2016 Landsat OLI imagery, and mean NDVI was calculated within each of the 15 sampling sites. Topographic variation was quantified as the standard deviation of altitude using data from the Shuttle Radar Topography Mission (SRTM) within each of sampling sites. Site extents for mean NDVI and topographic variation were defined by a 500-m buffer around each botanical plot noted above. Within-site conditions (υ) for each individual were assigned based on the site in which an individual was sampled, where all individuals sampled within a site were assigned the same within-site conditions. We initially evaluated nine 'within-site' models against a null model of isolation-bydistance (IBD) to evaluate the significance of within-site conditions alone in explaining current ruffed lemur gene flow. Within-site models included singular models for each of the seven within-site conditions, forest structure (a combination of canopy height, basal area, and stem density), and forest productivity (a combination of ENS, IVI, and NDVI). We combined within-site variables that performed better than IBD with resistance distance variables to evaluate the joint potential of within-and between-site metrics in explaining ruffed lemur gene flow. In total, we tested 15 composite models along with their component singular models. We log-transformed the response variable (D PS ; gene flow) and all predictor variables. Gravity models were performed using linear mixed effect models fit using Restricted Maximum Likelihood (REML), and model selection was conducted using log-likelihood (corroborated with Akaike Information Criterion; AIC). Finally, we used β estimates from all singular models as an evaluation of the variables' directional effect on connectivity, with positive β values suggesting within-site characteristics were attractive to ruffed lemur dispersal [99].

Resistance Analysis
We found a significant signature of isolation-by-Euclidean distance, which explained 4.4% (Pearson's r = 0.2091) of the observed population genetic structure ( Figure S1). Despite evidence of IBD, three surfaces-terrain ruggedness index (TRI), canopy height, and topographic position index (TPI)-were more strongly associated with genetic differentiation than geographic distance alone (Tables 1-3 and S1). Resistance to terrain was lowest in areas of low-to-moderate ruggedness where altitude did not vary greatly on a fine scale (below 10 ha; Figure 3 and Table 2). Relatively tall canopy (above 15 m) appears to better facilitate functional connectivity and gene flow than areas with low canopy height ( Figure 4 and Table 2). Furthermore, resistance was lowest on flat or moderate slopes leading to ridges ( Figure 5 and Table 2). The remaining five surfaces (1990 and 2016 forest cover, altitude, rivers, and fire density) explained slightly more variation than distance alone but were seldom chosen as the top model (less than 1.0% of the time) and thus were not considered in the composite analysis (Table 1). Although both 1990 and 2016 forest cover were rarely considered in the top model, results for both single surface resistance models indicated that ruffed lemur gene flow was more resistant through matrix than forest (Table 3). Additionally, resistance decreased monotonically with increasing altitude and exponentially with increasing fire density (Figures S2 and S3 and Table 2). Finally, rivers were predicted to cause more resistance to gene flow in ruffed lemurs than the intervening landscape (Table 3). Table 1. Results from bootstrap selection of optimized linear-mixed effects models on single surfaces. Rows shown in bold indicate models that performed better than Euclidean distance alone and were selected as the top model more than 1% (π ≥ 1.00) during 10,000 bootstrap iterations.         To assess the combined effects of terrain, canopy height, and topography, w four composite surfaces that combined TRI, canopy height, and TPI together in a tations and compared results against the isolated surfaces, as well as against str Euclidean distance. There was low correlation among all surfaces used in the c To assess the combined effects of terrain, canopy height, and topography, we created four composite surfaces that combined TRI, canopy height, and TPI together in all permu-tations and compared results against the isolated surfaces, as well as against straight-line Euclidean distance. There was low correlation among all surfaces used in the composite models, with Spearman rank correlation coefficients less than 0.10 in all cases (Table S2). There was no clear consensus on the top model predicting ruffed lemur gene flow, as both TRI and canopy height in singular, the composite of TRI and canopy height (Combination C), and the composite of TRI and TPI (Combination D) were all ranked similarly and were top models more than 10% of the time (Tables 4, 5 and S3). Resistance patterns in the composite surfaces showed an identical pattern for TRI as with the variable in singular, where resistance was lowest in areas of low-to-moderate ruggedness (Figures 6 and 7). Canopy height was transformed to distance when combined with TRI, meaning this composite surface was essentially identical to TRI in isolation (Combination C; Figure 6). However, for the other two composite models containing canopy height (Combinations A and B), we found the same pattern as canopy height in singular, where taller trees above 15 m appear to facilitate dispersal (Figures S4 and S5). Transformations for TPI showed greater variability, although for composite surfaces containing TPI resistance was lowest on flat or moderate slopes leading to ridges, as was found in the singular surface (Figure 7, Figures S4 and S5). For both highly ranked composite surfaces (Combinations C and D), TRI contributed to more than 95% of the surface, and results were therefore driven primarily by TRI (Table S4). This suggests that TRI on its own was the top model approximately 60% of the time (a sum of TRI in singular and contributions from Combinations C and D), lending support for TRI as the best predictor of the genetic data. However, canopy height was also a top-ranked model nearly 20% of the time and, therefore, cannot be excluded as a strong predictor (Table 4). Our visualization of the optimized TRI model shows a high degree of current (gene) flow throughout the eastern extent of our study area, particularly in the eastern portion of the remaining forest ( Figure 8). Our visualization of the canopy height model revealed a relatively uniform current flow throughout much of the remaining rainforest in our study extent, although the visualization does suggest current flow may be greatest in the eastern and western boundaries of the remaining forest ( Figure 9). Together, these results suggest a high degree of connectivity throughout all localities sampled, although the northern-most site (VIAVY) seems to have weaker current flow than any of the other sites (Figures 8 and 9).    For both highly ranked composite surfaces (Combinations C and D), TRI contributed to more than 95% of the surface, and results were therefore driven primarily by TRI (Table  S4). This suggests that TRI on its own was the top model approximately 60% of the time (a sum of TRI in singular and contributions from Combinations C and D), lending support for TRI as the best predictor of the genetic data. However, canopy height was also a topranked model nearly 20% of the time and, therefore, cannot be excluded as a strong predictor (Table 4). Our visualization of the optimized TRI model shows a high degree of

Gravity Model
When investigating within-site drivers of dispersal, we found equal support for the null model of IBD and the test models of IVI and NDVI, with no support for the other test models (Table 6). Unsurprisingly, we found a negative relationship between distance and D PS , suggesting fewer shared alleles between individuals at increasing geographic distances (β = −0.030, SE = −0.003, df = df = 3685, t = −9.36, p < 0.001). We found a significant and slightly positive relationship between the average IVI of Varecia food resources and D PS , suggesting that these food resources were attractive to ruffed lemur dispersers (β = 0.093, SE = 0.036, df = 60, t = 2.62, p = 0.011). Finally, we found a negative (but non-significant) relationship between mean site NDVI and D PS , suggesting that sites with higher productivity values have a reduced attractiveness to dispersal than sites with lower NDVI values (β = −0.261, SE = 0.465, df = 60, t = −0.56, p = 0.576). Together, these results indicate that sites such as TALA, SAKA, and HG are more likely to attract dispersers due to their relatively higher values of average IVI and lower values of average NDVI compared to sites such as MGV, MAND, and MALA, which display lower average IVI and higher average NDVI ( Figure 1 and Table S5). Finally, we combined the best-supported within-site models (IVI and NDVI) with the two best-supported resistance models (TRI and DSM) to evaluate the combined influence of within-and between-site drivers on ruffed lemur dispersal. First, all models containing one or both resistance variables were more strongly supported than the null model of IBD (Table 7). Two models-canopy height in singular and canopy height combined with NDVI-showed equal support as the best models predicting dispersal in ruffed lemurs throughout the Ranomafana region ( Table 7). As expected, in the composite model, we found a significant and negative relationship between canopy height resistance distance and D PS , confirming that greater resistance distances reduced gene flow (β = −0.216, SE = 0.024, df = 3683, t = −9.01, p < 0.001). Furthermore, we found a negative (but non-significant) relationship between mean site NDVI and D PS , similar to the variable in singular (β = −0.488, SE = 0.437, df = 60, t = −1.12, p = 0.269). As with the composite model, we found a significant and negative relationship between canopy height resistance distance and D PS when evaluating the variable in singular (β = −0.215, SE = 0.024, df = 3683, t = −8.96, p < 0.001). Table 7. Results for composite, within-site, and between-site gravity model predictions. Rows shown in bold indicate the best supported model based on log-likelihood criteria and corroborated using AICc.

Discussion
The main objective of this study was to evaluate the between-and within-site environmental drivers of black-and-white ruffed lemur dispersal throughout Madagascar's southeastern Ranomafana region. We found that the predominant between-site drivers of ruffed lemur dispersal were terrain ruggedness index (TRI) and canopy height, where areas of low-to-moderate ruggedness and relatively tall canopy heights (above 15 m) facilitated functional connectivity and gene flow. The other six environment variables evaluated (topographic position index, 1990 and 2016 forest cover, altitude, rivers, and fire density) were not significant predictors of ruffed lemur dispersal, although resistance was lowest within forests, at moderate slopes, at lower altitudes, in areas with lower fire density, and outside of rivers. Productivity metrics, including the Normalized Difference Vegetation Index (NDVI) and the average Importance Value Index (IVI) of ruffed lemur food trees, were the most influential within-site drivers of ruffed lemur dispersal, with animals being more likely to disperse into areas with lower NDVI and higher IVI values. However, when between-and within-site drivers were combined, models containing only between-site resistances were consistently the best supported. We, therefore, conclude that the within-site environmental variables tested were less influential in ruffed lemur dispersal decisions than between-site factors.

Between-Site Influences on Ruffed Lemur Dispersal
The best-supported between-site drivers of ruffed lemur dispersal in the Ranomafana region were terrain ruggedness index (TRI) and canopy height. Ruffed lemur dispersal was facilitated by areas of low-to-moderate ruggedness and impeded by increasingly rugged terrain. Our results add to a growing body of evidence identifying areas of high ruggedness as impediments to gene flow (including lizards, Sceloporus occidentalis: [100]; stone marten, Martes foina: [101]; Northern quoll, Dasyurus hallucatus: [102]; and tigers, Panthera tigris: [31]; but see Balkenhol et al. [49], wolverine, Gulo gulo). For instance, tigers preferentially selected the least rugged dispersal paths from within their preferred landscapes-forested areas with low anthropogenic presence-suggesting that avoiding areas of high terrain ruggedness was secondary to moving through their preferred landscape [31]. We found a similar pattern in ruffed lemurs, with gene flow predicted to be highest in areas of low-tomoderate ruggedness in the eastern boundary of the remaining forest of the Ranomafana region (Figures 3 and 8). Interestingly, our results predicted that ruffed lemurs would avoid areas of the lowest ruggedness ( Figure 3). The reasoning for this may be twofold: low ruggedness areas of the eastern Ranomafana region are associated with increased anthropogenic presence and are preferentially used by humans for both irrigated and shifting agriculture (i.e., flooded rice paddies and slash-and-burn or tavy; [103,104]), whereas low ruggedness areas in the western region are associated with higher altitudes and shorter forest canopies [65][66][67]. In this way, similar to tigers, the ruffed lemurs in our study area are likely dispersing through the least rugged terrain that enables them to avoid anthropogenic landscapes while also remaining within their preferred tall-canopied habitats [31]. Alternatively, the relationship found between topography and dispersal may be driven by ruffed lemur behavioral ecology, similar to what has been hypothesized in wolverines [49]: ruffed lemurs segment their territories along ridgelines (Baden, pers. comm.), making it possible that individuals preferentially move along ridgetops with less rugged terrain to avoid core territories of unknown individuals.
In addition to terrain ruggedness, canopy height was also predicted as a significant driver of ruffed lemur dispersal in the Ranomafana region. We found that relatively tall canopy (above 15 m) facilitated and short canopy height impeded dispersal between localities. In other forest-dwelling species, canopy height or structural homogeneity were also found to be significant predictors of gene flow (Western capercaillie, Tetrao urogallus: [105]; black-eared mouse, Peromyscus melanotis: [106]), although unsurprisingly the forest-dwelling scarlet macaws (Ara macao)-which disperse via flight-were not impacted by either [107]. Ruffed lemurs rely on tall, broad-canopied fruiting trees due to their obligate frugivory [53,56] and, despite their ecological flexibility, are found at greater densities in localities characterized by these features [57,108,109]. It is thus unsurprising that taller forest canopy facilitates dispersal, as animals tend to preferentially move through landscapes that are most similar to their preferred habitat [11,12]. The high rates of gene flow identified by our study likely correspond with areas of tall canopy at lower altitudes along the eastern boundary of the remaining forest in the Ranomafana region, similar to patterns found in other parts of the world [65][66][67] (Figure 9). Our results also predicted high rates of gene flow throughout fragmented forests along the western edge of the region ( Figure 9); however, dispersal through these areas is unlikely given its higher altitude, as there are few records of ruffed lemurs occurring above 1200 m and no evidence above approximately 1350 m [110][111][112].
Curiously, forest cover was not identified as a significant driver of ruffed lemur dispersal, as in Baden et al. [46], which may reflect scale-dependent differences in dispersal drivers. Specifically, the simple presence of forest may be an important driver of gene flow for longrange and multigenerational dispersal, while the quality of forest may be most influential for short-range dispersal decisions. Alternatively, these results may simply mirror those found by Milanesi et al. [105], where forest structure derived from LiDAR remote-sensing yielded better estimates of gene flow compared to traditional land cover data.

Within-Site Influences on Ruffed Lemur Dispersal
Ultimately, within-site habitat features in this study were less influential than betweensite features in driving ruffed lemur dispersal in the Ranomafana region. These findings are similar to other recent studies using similar methodologies (black-capped vireos, Vireo atricapilla: [112]; Arizona treefrog, Dryophytes (Hyla) wrightorum: [113]). Although not significant, productivity measures-Normalized Difference Vegetation Index (NDVI) and the average Importance Value Index (IVI) of ruffed lemur foods-were the best supported habitat features driving ruffed lemur emigrations and immigrations. Sites with relatively lower IVI and higher NDVI, such as Mangevo (MGV), Mandiandry (MAND), and Malazamasina (MALA), likely act as source populations (i.e., source of emigrants), and those with relatively high IVI and lower NDVI, such as Talatakely (TALA), Sakaroa (SAKA), and Harangangavo (HG), are likely most attractive to dispersers and may act as sinks (Figure 1; Table S5). This finding is further supported by the presence of several first-generation immigrants identified in both TALA and SAKA, suggesting that these sites are indeed attractive to dispersers [59].
There is high within-site structural variability between the sampling localities evaluated in this study, with more pristine and higher quality sites such as Mangevo (MGV), Tandrokaomby (TAND), and Amboasary (AMBO) having taller canopies, larger basal area, and greater stem density than disturbed, lower quality sites of Talatakely (TALA) and Sakaroa (SAKA ; Table S6) [59]. Many forest-dwelling species prefer to move through higher quality, complex forest stands, as opposed to open habitat [114][115][116], and preferentially incorporate these stands into their home range [116][117][118]. Furthermore, several studies have found that higher quality sites (relative to the species' biology) facilitated gene flow and may act as sources of emigrants for the species (American martens, Martes americana: [52]; Columbia spotted frogs, Rana luteiventris: [98]; Blainville's horned lizard, Phrynosoma blainvillii: [119]). It was, therefore, surprising that within-site forest structure was not identified as a significant predictor of ruffed lemur dispersal. It is possible that our metrics did not fully capture details of forest structure that are important to ruffed lemur movement. In the future, LiDAR or photogrammetry could enable a more detailed quantification of forest structure that our study may have missed, while allowing researchers to measure other important characteristics such as canopy density and size and distribution of canopy gaps [115,[119][120][121].

Spatial Variation in Response to Environmental Variables
In addition to our evaluation of between-and within-site drivers of dispersal, we compare our regional results to a recent range-wide assessment by Baden et al. [46] to evaluate how ruffed lemur dispersal might vary across spatial scales. The study by Baden et al. [46] assessed pattens of long-range dispersal and multigenerational gene flow and found that long-range ruffed lemur dispersal is facilitated primarily by available forest cover. By contrast, the present study evaluated smaller-scale, more regional gene flow from typical dispersal events and found weak evidence for this same relationship. Forest cover is contiguous throughout our study area, making it possible for ruffed lemurs to avoid the matrix when dispersing through the region, perhaps explaining why 1990 and 2016 forest do not appear to play significant roles in regional ruffed lemur dispersal. We did, however, find the height of contemporary forest-a strong correlate of 2016 (contemporary) forest cover (ρ = 0.77; Table 6) and a moderate correlate of 1990 (historic) forest cover (ρ = 0.44; Table S2)-was a significant predictor of dispersal, with taller canopy facilitating dispersal. Tall canopy in the Ranomafana area is often related to higher quality forest [60,65] where ruffed lemurs are found at their highest abundances [53,108,109]. Additionally, ruffed lemurs rely on large, broad canopy trees for quadrupedal movement, access to sufficient fruit [53], and successful reproduction [122]. Therefore, ruffed lemurs' utilization of tall canopy trees for dispersal is expected, particularly as animals preferentially use areas most similar to their chosen habitat for dispersal [11,12].
Proximity to human settlements was found to be the main deterrent to long-range ruffed lemur dispersal and multigenerational gene flow [46]. In comparison, we did not find a significant influence of fire density-a more nuanced metric of anthropogenic presence-on short-range ruffed lemur dispersal. Although not significant, we did find an exponential increase in resistance with increasing fire density, similar to the pattern found by Baden et al. [46] with proximity to human settlements. Topography was not found to be a significant driver of long-range ruffed lemur dispersal across their remaining range [46], although we did find the terrain ruggedness index to be a significant betweensite driver of short-range dispersal in the Ranomafana region, with low-to-moderate ruggedness facilitating dispersal. Finally, neither study found any influence of rivers on dispersal, despite the role rivers have played in the biogeography of Malagasy primates [110,123]. The permeability of rivers is likely dependent on several factors, including size, flow, and access to bridging structures (e.g., overhanging trees or anthropogenic structures such as bridges), with larger waterways expected to be less permeable than minor rivers [124,125]. Only minor waterways are present within our study and, therefore, unlikely to pose dispersal barriers to the species, as they are likely able to disperse across these minor rivers via natural canopy bridges.

Conclusions
In summary, we found that between-site environmental characteristics-specifically, terrain ruggedness and forest canopy height-were better predictors of short-range ruffed lemur dispersal than within-site characteristics such as NDVI and Importance Value Index of ruffed lemurs' foods in the Ranomafana region. Furthermore, we found evidence of scale-dependent environmental influences on ruffed lemur dispersal, with topography only driving short-range dispersal decisions. Forest characteristics influenced dispersal at all scales, with forest presence facilitating multigenerational, long-range dispersal and taller forest heights facilitating short-range dispersal in ruffed lemurs. Together, results from this and earlier studies [46] highlight the importance of high-quality forest to sustaining ruffed lemur gene flow across spatial scales. Given the accelerating forest modification, loss, and fragmentation throughout Madagascar's eastern rainforest escarpment [55,112,[126][127][128][129], it is likely ruffed lemur dispersal, along with the dispersal patterns of countless other arboreal taxa [130], will become increasingly restricted in the future. Dispersal capacity and limitations are primary drivers of community structuring across primate taxa at regional spatial scales, more so than niche differentiation [3,130]. If ruffed lemur dispersal becomes significantly limited by forest loss, this may lead to more divergent primate community structure, even between neighboring localities [130]. Shifting community structure, decreasing forest suitability, and limited dispersal ability will have cascading effects, and will likely result in massive decreases in population sizes beyond what would be expected with shifting suitable habitat alone [131,132]. Therefore, retaining high-quality forests (particularly in areas without strict protection) and forest connectivity is paramount to retaining gene flow of ruffed lemurs within their southern range, ultimately buffering against ongoing population declines and local extinctions in this critically endangered species.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/genes14030746/s1, Methods S1: Forest Cover Classifications [126,[133][134][135][136]; Figure S1: Linear regression between pairwise Euclidean and genetic (Ar) distances; Figure S2: Resistance transformation plot for single surface altitude; Figure S3: Resistance transformation plot for single surface fire density; Figure S4: Resistance transformation plot for 'Composite A' surface of (A) Terrain Ruggedness Index (TRI), (B) canopy height, and (C) Topographic Position Index (TPI); Figure S5: Resistance transformation plot for 'Composite B' surface of (A) Terrain Ruggedness Index (TRI) and (B) canopy height; Figure S6: Two-dimensional visualizations of unmixed feature space for High Albedo, Low Albedo, and Vegetation end members from 2009 mosaic; Figure S7: Customized decision tree for thematic classifications. Bands used for classification included b1: high albedo end member; b2: vegetation end member; and b3: low albedo end member; Table S1: Parameter estimates from mixed effects model fit to the optimized singular resistance surfaces; Table S2: Spearman rank correlation between singular resistance surfaces; Table S3: Parameter estimates from mixed effects model fit to the optimized composite resistance surfaces; Table S4: Percent contribution of each for each singular surface comprising the optimized composite surfaces; Table S5: Mean height, basal area (m 2 /ha), stem density (stems/ha), effective number of species (ENS; exp(H')), mean IVI of Varecia food species, mean NDVI, and standard deviation of altitude at 15 sites in the Ranomafana region; Table S6: Landsat scene information from composite scenes used in generating thematic maps.