Effects of Forest Fragmentation on Connectivity and Genetic Diversity in an Endemic and an Invasive Rodent in Northwestern Madagascar

1. Introduction

The island of Madagascar is a globally important biodiversity hotspot with high rates of endemism [1,2,3]. Yet habitat loss, mainly via deforestation, poses an urgent threat to the conservation of biodiversity, in particular to endemic species on the island [4,5]. Since the 1950s, Madagascar has lost approximately 44% of its total forest area leading to a patchwork of habitat fragments exposed to edge effects and surrounded by non-habitat matrix [6]. Generally, as habitat area decreases, community-level species richness and the abundance of individual species will decrease [7,8,9]. Following habitat loss in fragmented landscapes, a greater proportion of remaining habitat is subject to edge effects which can threaten populations of edge-sensitive species, a potentially confounding factor as habitat loss and increasing edge ratios occur synergistically [10,11,12]. Alternatively, edge-tolerant or ecologically flexible species may see positive outcomes following processes such as forest loss as more diverse habitat types become available [10,13]. Malagasy mammals do not respond to forest loss and fragmentation in uniform ways [14,15]. Species-specific responses to edges [16,17], vagility between remaining habitat patches [18,19] and competition with invasive species [20,21] are important research topics relevant to researchers and conservation practitioners alike in Madagascar.

When organisms are subjected to habitat loss and fragmentation, genetic and genomic methods are useful to explore and quantify negative outcomes on population genetic composition [22,23]. Populations in fragmented landscapes may be subject to decreased gene flow between isolated subpopulations which can lead to decreased genetic diversity across a landscape and increased genetic drift in isolated fragments [24]. These processes can lead to loss of genetic diversity through drift and inbreeding and eventually to the fixation of deleterious alleles within populations [25,26]. Increasing proportions of edge habitats in fragmented landscapes may further affect the connectivity and genetic diversity of organisms. For example, Radespiel et al. [27] found local genetic differentiation in a small-bodied lemur species, Microcebus ravelobensis, between edge and interior habitats within the same forest. However, the negative effects of smaller and more isolated populations may be mitigated if connectivity is high within landscapes and organisms are able to disperse across matrix and edge zones between habitat patches [28,29]. After processes such as deforestation, the new habitat types that arise in the matrix may facilitate dispersal of some organisms, in particular ecologically flexible and/or invasive species [30].

To evaluate the effects of habitat fragmentation, our study focused on two species of rodents in northwestern Madagascar: the endemic western tuft-tailed rat (Eliurus myoxinus), and the invasive black rat (Rattus rattus). As a geographically widespread species, E. myoxinus are found throughout the western and northern dry deciduous forests [31,32]. They are a small-bodied (~60 g), solitary, nocturnal mammal with a primarily arboreal lifestyle, use nests, burrows, and tree holes as sleeping sites, and consume mainly fruits [33,34,35]. Other basic natural history information such as dispersal patterns and distances are still lacking for E. myoxinus. As with other Malagasy fauna, studies on forest loss and fragmentation have found mixed or unclear effects. Species-specific and regional differences appear to play a large role in determining the direction and scale of effect, rather than a single unifying ecological pattern. For example, capture rates of Eliurus spp. increased with increasing forest area in some dry forest sites [20,33]. In other dry forest sites no evidence for an area effect was found but signals of decreased abundance in fragmented forests were observed in relation to increasing distance from habitat edge [21]. Conversely, populations of Eliurus spp. in the North and West of Madagascar display relatively consistent biogeographic patterns, with weak signals of genetic differentiation from large-scale geographic barriers such as rivers, and a lack of strong isolation by distance pattern suggesting that ecological factors such as habitat preference and elevational gradients play a more important role in the biogeography of these rodents [32,36,37,38].

In comparison, R. rattus is a larger (~100 g), more gregarious species with a generalized diet that is found in a variety of habitats from urban and rural human settlements to primary and secondary forests [39]. When enough resources are present, R. rattus do not disperse over long distances, preferring to stay in close proximity to kin [40]. The introduction of R. rattus in Madagascar very likely occurred recently, within the last 1,000 years, in port cities like Tôlanaro in the east or Antsiranana in the north via trading routes from Africa and Asia [41,42]. A complicated picture of fragmentation effects has so far been documented in R. rattus, perhaps due to the complicated nature of defining habitat preference for an invasive generalist species. Capture rates of R. rattus declined with increasing forest area in western dry forest sites [20], while forest area did not affect capture rates in other regions in the northwest and southwest [21,33]. In Central Madagascar, patterns of genetic differentiation in R. rattus were only weakly affected by habitat heterogeneity, and significant patterns of isolation by distance were detected as is expected in an invasive generalist [43]. However, Brouat et al. [42] found an unexpectedly weak signal of geographic distance for R. rattus populations, and suggested dispersal alongside humans in a non-euclidean pattern, such as via boats along coastlines, was responsible for the spatial genetic patterns observed. Several studies have investigated the potential for competition between introduced and endemic rodent species, but these studies found no effect or inconclusive results [21,33]. Even if direct competition is not occurring between invasive and endemic rodents, elucidating what ecological factors allow for each of these species to persist in fragmented forest landscapes constitutes an important question from an ecological and conservation perspective.

The goal of our study was to determine the molecular effects of habitat loss and fragmentation in these two contrasting species of rodents, E. myoxinus and R. rattus, in two fragmented forest landscapes in northwestern Madagascar. Given one arboreal, relatively specialized, and forest-dependent endemic species and one ecologically flexible species with a generalized diet that has adapted to a variety of habitat types, different outcomes can be expected for these two species. We aimed to answer three main research questions for each of these two species in order to compare and contrast their movement patterns and the resulting population connectivity when living in fragmented landscapes: 1) Do endemic and invasive rodents use non-forest matrices for local movement or dispersal? 2) Are connectivity and genetic structure within these landscapes affected primarily by isolation by distance or by habitat-related factors? 3) Which ecological factors affect the genetic diversity of subpopulations living in fragmented landscapes?

2. Materials and Methods

2.1. Study Site and Sampling

We systematically captured animals across two landscapes in northwestern Madagascar (Mariano Classified Forest: MCF; and Ankarafantsika National Park: ANP) that contain a network of forest patches ranging in area from 0.8 to 114.6 ha with adjacent continuous forest sites ranging from 3,683 to 130,390 ha (Figure 1 and Table S1). Both landscapes are part of the dry deciduous forest ecoregion in Madagascar, an area of global conservation concern [44]. The two landscapes are separated by approximately 100 km of anthropogenically disturbed area containing settlements, rice agriculture, and grassland. Sampling occurred between June and October 2017 and 2018 during the dry season in Madagascar. While broadly similar in terms of habitat type and climate, the studied landscapes differed in several factors. The spatial scale at which we studied each landscape was different; edge-to-edge geographic distance between fragments varied from 30 m to 18,182 m in MCF and from 20 m to 6,342 m in ANP. In MCF, the forests were situated at low altitudes (10 - 100 m above sea level, a.s.l) and were more humid than in ANP due to higher amounts of surface water and proximity to the ocean. In ANP (90 - 208 m a.s.l), no surface water was available during the dry season, and fragments were surrounded by a relatively homogeneous matrix of grassland savannah which has been present in its current configuration for at least 70 years [19]. In contrast, fragments in MCF were surrounded by a heterogenous matrix containing varying densities of palm trees (Bismarckia nobilis) and a generally anthropogenic landscape featuring villages, rice fields, and roads [45]. ANP is managed by Madagascar National Parks (MNP), while MCF is managed by the Vondron’Olona Ifotony (VOI), a community-based forest protection program.

In ANP, we established 31 transects in 27 forest fragments and four transects in continuous forest sites. In MCF, we established 18 transects in 13 forest fragments and eight transects across three patches of continuous forest (Figure 1). All transects were established perpendicular to the forest edge using a compass and a handheld GPS unit during installation. Forest edges and surface areas were defined according to Steffens and Lehman [19]. Polygons for each fragment were created by walking the perimeter of the fragment with a handheld GPS. Transects ranged in length from 40 m to 490 m into forest fragments (the shortest distance to reach the center of the fragment) and up to 1,000 m in the continuous forest. We set Sherman live traps (Sherman Traps, Inc., Tallahassee, FL, USA) baited with bananas in pairs every ten meters at approximately one meter above the ground along these transects. Captured animals were anesthetized and then marked (via ear cuts), weighed, sexed, and a small biopsy was taken from the ear, placed in Queen’s Lysis Buffer, and stored in a fridge. For detailed description of trapping methodology and results see Andriatsitohaina et al. [21] and Kiene et al. [46] for anesthetic procedures. We obtained an additional six samples of E. myoxinus collected in two other sites (Andoharano, Ambanjakely) in ANP in 2004 for comparative spatial analyses, see Rakotondravony and Radespiel [47] for further details on these study sites. We used QGIS (version 2.18.9) to determine the distance to the next edge, distance to the next continuous forest, and forest area for each individual animal captured.

All field protocols were reviewed and approved by the appropriate ethics committees at the University of Toronto, Canada (Protocol # 20011950) and the Institute of Zoology, University of Veterinary Medicine Hannover, Foundation, Germany. All fieldwork adhered to the legal requirements of Madagascar and were approved by the Ministère de l’Environnement, de l’Ecologie et des Forêts and Madagascar National Parks (Permits # 81/17/MEEF/SG/DGF/DSAP/SCB.Re, #151/17/MEEF/SG/DGF/DSAP/SCB.Re, #82/18/MEEF/SG/DGF/DSAP/SCB.Re, #92/18/MEEF/SG/DGF/DSAP/SCB.Re).

2.2. Labwork and Data Preparation

We used restriction site-associated DNA sequencing (RADseq) to generate genomic data for the population genomics analyses [48]. RADseq generates sequence datasets for a large number of sites across the genome and is generally useful for the study of non-model organisms as no prior genetic information is needed [49]. DNA was extracted from ear biopsies using a DNeasy Blood & Tissue Kit (Qiagen) following an optimized version of the manufacturer’s protocol. DNA concentrations were estimated using the Fluorometer dsDNA HS Assay Kit (Qubit) to select samples of a sufficiently high concentration. We obtained 93 samples of E. myoxinus (n = 77 individuals from ANP, n = 16 individuals from MCF) and 106 samples of R. rattus (n = 49 individuals from ANP, n = 57 individuals from MCF) for RADseq.

RADseq libraries were prepared using 100-200 ng of genomic DNA and the TruSeqNano DNA HT kit (Illumina). The restriction enzyme SbfI was used to generate DNA fragments. Libraries were prepared in sets of 24 samples sorted based on their measured DNA concentration. Each sample was ligated to one of 48 available P1 adapters with a unique 5 base pair molecular identifier allowing individual identification of the samples after sequencing. All fragments were randomly sheared resulting in fragments with an average size of 550 bp. We ligated sheared DNA fragments to the P2 adapter, and all fragments with both P1 and P2 adapters were amplified in 10 polymerase chain reaction cycles. We verified DNA concentration and fragment sizes of the amplified libraries using a Qubit Fluorometer and a Bioanalyzer 2100. Finally, we sequenced the libraries using 150bp paired-end reads on an Illumina HiSeq3000 at the GeT-PlaGe platform (Toulouse, France).

Due to a lack of a reference genome for E. myoxinus, we used a de novo alignment method in Stacks 2.5 [50] to call SNPs. We first filtered out low quality reads using process _radtags, clone_filter and kmer_filter. We created a catalog and aligned our samples using the denovo _map pipeline in Stacks utilizing ustacks, sstacks, and tsv2bam. We set m, M, and n parameters to two, seven, and eight after optimization using protocols from Rochette and Catchen [51]. Lastly, we ran gstacks and populations retaining only loci found in at least 80% of samples (r = 0.80), with a minor allele frequency of at least 5% (min_maf = 0.05), a maximum heterozygosity of 0.7 (max_het = 0.7), and we retained only one random single nucleotide polymorphism (SNP) per RAD locus (--write_random_snp) in order to avoid linkage disequilibrium in downstream analyses. After filtering we retained 66 samples of E. myoxinus, 52 from ANP and 14 from MCF (Table S2).

We used a similar protocol and parameters for R. rattus whenever possible, but with a reference genome available. In Stacks 2.5 [50], we used process _radtags, clone _filter, and kmer _filter to filter out low-quality reads. We then aligned our samples of R. rattus to the reference genome of congener R. norvegicus (Rnor_6.0, GCA_000001895.4) using the Burrows-Wheeler Alignment (BWA) with BWA-MEM [52]. We then ran gstacks and removed read pairs of the same length (--rm-pcr-duplicates). We ran populations retaining only loci found in at least 80% of samples (r = 0.80), with a minor allele frequency of at least 5% (min_maf = 0.05), a maximum heterozygosity of 0.7 (max_het = 0.7), and we retained one random single nucleotide polymorphism (SNP) per RADtag (--write_random_snp). After filtering, we retained 81 samples of R. rattus, 43 from ANP and 38 from MCF (Table S3).

As a first test to assess if each the two landscapes contained distinct populations, we ran ADMIXTURE [52]. We tested a varying number of ancestral populations (K=2-10) and calculated cross-validation errors (CV) for each K-value to determine the most likely K. We then used the CLUMPAK program to infer cluster assignments across 20 independent runs of ADMIXTURE [54] which revealed that each landscape contained a distinct population with comparably low levels of admixture between them (E. myoxinus: < 30%, R. rattus: < 27%) (Figure 1, Table S3). Therefore, all subsequent analyses were performed separately for each species in each landscape. Due to limited sample sizes once species were separated at each landscape, no additional sub-analyses were completed (such as by sex or fragmented and continuous populations).

2.3. Analyses of Local Movement and Dispersal

We calculated a coefficient of genetic relatedness for each dyad of individuals for each landscape separately to investigate the distribution of related individuals in space. Distribution of related individuals reflects well the dispersal patterns within a species and can be used to detect potential dispersal events [55]. We generated pairwise relatedness by multiplying the kinship coefficient from the relatedness2 function in VCFtools [56] by two to assign dyads into kinship classes based on these r values. Classes were defined as first-order relatives (r > 0.375), second-order relatives (0.375 > r > 0.1875), or distantly related individuals (0.1875 > r > 0) in congruence with Combs et al. [57]. Dyads with r ≤ 0 were considered unrelated. We calculated geographic distances between all related dyads and calculated the mean, median, and range of geographic distances for each kinship class. We ran Mann Whitney U tests in R to test for significant differences in geographic distances between second-order relatives and distantly related dyads within and between each species. First-order dyads were not statistically analysed in both species due to low sample sizes.

We additionally created a genetic dissimilarity matrix using the program BED2DIFFS_V1, a method that calculates average genetic dissimilarity based on allelic differences observed between samples [58] and generated a geographic distance matrix (Ersts, American Museum of Natural History, Center for Biodiversity and Conservation, http://biodiversityinformatics.amnh.org/open_source/gdmg). We used these two matrices to create Mantel correlograms using the mantel.correlog function in the Vegan R package with 999 permutations [59]. Mantel correlograms test for spatial autocorrelation between genetic and geographic distances in a predefined set of distance classes. We selected distance classes of 200 m width due to the expectation for short dispersal distances in R. rattus [60] and used the same distance classes for E. myoxinus for comparability. Only dyads from individuals of less than half the maximum Euclidean geographic distance were included in this analysis to reduce a potential bias due to individuals from the edges of the sampling space [61].

2.4. Analyses of Connectivity and Genetic Structure on the Landscape Level

In order to assess the effect of isolation by distance (IBD) across the entire population in each landscape we first conducted Mantel tests using the genetic dissimilarity matrix and the geographic distance matrix described above. We ran Mantel tests with 999 permutations on each species and landscape separately using the Vegan R package [59].

We used three complementary methods to determine sub-population connectivity and genetic structure within the two landscapes. For all three analyses we used a slightly reduced dataset with one individual of each closely related dyad being removed (n = 4). First, we ran ADMIXTURE [53] to assess genetic structure via ancestry estimation across each landscape using the same methods described above on each species and landscape independently. We pooled population ancestry estimates from ADMIXTURE for all individuals from each forest site into one pie chart to visualize population structure across the landscape.

Secondly, we used the Estimating Effective Migration Surfaces (EEMS) software package to investigate migration rates and potential barriers to gene flow across each landscape. EEMS uses a stepping-stone model to estimate effective migration rates between adjacent subpopulations (demes) and produces a visual migration surface of estimated genetic barriers and connectivity [58]. Specifically, the program highlights geographic areas where genetic distance between samples increases more quickly (i.e., low genetic connectivity) or more slowly than expected (i.e., high genetic connectivity). Following Barratt et al. [62] we used the genetic dissimilarity matrix generated by BED2DIFFS_V1 and raw geographic coordinates for each sample as input files for EEMS. We ran EEMS using deme sizes of 1,000 (the spatial resolution of the analysis) after evaluating the effect of varying resolution (250, 500, and 750, results not shown). Running EEMS with various deme sizes showed similar relationships across the landscapes and species, but due to their highest spatial resolution, we only present results from the runs with 1,000 demes. We produced migration surfaces for each species on each landscape separately and produced additional migration surfaces for R. rattus for the Northern and Southern patches of forest in MCF due to the large spatial scale of the site and sufficient sample size. Our analyses used two MCMC chains, each with a length of 25,000,000 iterations, with a burn-in of 1,000,000 iterations. We examined the log posterior between the chains and combined results from both using the rEEMSplot package. We then plotted surfaces of effective migration rates across each landscape to identify geographic barriers and areas of connectivity.

Thirdly, we ran the ResDisMapper program [63], which estimates resistance to dispersal across landscapes. ResDisMapper calculates IBD residuals for each dyad of individuals and maps resistance values based on their deviations from the general IBD trend. Genomic data were input into ResDisMapper using the GENEPOP format, and a Euclidean distance matrix was used to produce IBD residuals using the “relative similarity method” (Method 1). We ran ResDisMapper with default values of a final resolution of 50, confidence intervals set to 0.95, for 1,000 repetitions. The significance threshold for resistance contours was set to 0.05. For MCF we set the maximum distance for dyads at 10,000 m to ensure a continuous distribution of IBD values [63] as sampling did not occur continuously across MCF. Our samples from ANP only included dyads less than 5,000 m apart sampled continuously across the landscape. We used the rdm_mapper function to produce a resistance map which was then overlaid on a map of each landscape produced in QGIS.

2.5. Ecological Determinants of Genetic Diversity

We calculated three parameters of genetic diversity to determine which ecological factors affected genetic diversity of the two rodents across the landscapes. The three individual-based parameters were the inbreeding coefficient (F_IS), observed heterozygosity (H_O), and expected heterozygosity (H_E). These three parameters were calculated using the VCFtools --het command [56] for each individual E. myoxinus and R. rattus. Generalized additive models (GAMs) were built for each of these dependent variables (Link function = identity, family = Gaussian) using maximum likelihood estimation in the R package mgcv [64] and the fixed factors species (E. myoxinus, R. rattus), region (ANP, MCF), sex (male, female), forest type (continuous, forest fragment), body mass (a proxy for body condition), forest area (of the fragment where an individual was captured), distance to edge (of an individual to the nearest forest edge), distance to continuous forest (of a fragment edge to the nearest continuous forest edge, a proxy for isolation) and the percentage of edge habitat (area of the fragment within 50 m of the forest edge). Normality of all continuous variables was assessed using Shapiro Wilk’s tests [65]. If a variable was not normally distributed, it was transformed using either log or square root transformations and Q-Q plots were visually inspected to ensure improvement. We tested for correlations between all continuous predictor variables (Body Mass, Forest Area, Distance to Edge, Edge Percentage, and Distance to Continuous Forest) to ensure independence of the fixed factors included in a given model (Table S4). We modeled each genetic diversity parameter (a) with the complete dataset (n = 123 individuals) and (b) just with the individuals caught in forest fragments (n = 90 individuals) to also test for variables that are not relevant for continuous forest sites (i.e., ‘percentage edge habitat’ and ‘distance to continuous forest’). The separation of correlated variables into different models resulted in two GAMs built for the overall data set (Models A-B) and three GAMs built for the fragment-only data set (Models C-E, Table 1). Additionally, we ran each model for each species separately to control for the potentially confounding interaction of species with other variables.

We determined the best model for each full model A-E using a model-averaging algorithm implemented in the “dredge” command from the R package MuMIn [66]. This process explores all possible combinations of factors in a given model and ranks them according to Akaike’s Information Criterion value [67]. We considered as best models those with a low corrected Akaike’s Information Criterion (ΔAICc < 2) and the highest AICc weight (ω) [67]. We used local regression smoothing and Loess smoothing to visualize nonlinearities in the models [68].

3. Results

3.1. Genomic Dataset

After filtering with the STACKS pipeline, we extracted 2,863,702 loci (42,855 – 59,901 per individual) with a mean depth of 9.79x coverage (5.05x – 23.82x per individual) for E. myoxinus individuals (n = 52) from ANP, and 917,682 loci (53,070 – 69,657 per individual) with a mean depth of 10.09x coverage (8.57x – 12.39x per individual) for E. myoxinus individuals (n = 14) from MCF. A total of 2,005,917 loci (40,542 – 49,097 per individual) with a mean depth of 9.68x coverage (5.40x – 19.43x per individual) were recovered for R. rattus individuals (n = 43) from ANP, and 1,761,411 loci (34,451 – 48,855 per individual) with a mean depth of 10.42x coverage (5.90x – 15.04x per individual) were found for R. rattus individuals (n = 38) from MCF (Table S2).

3.2. Patterns of Local Movement and Dispersal

We identified 26 related E. myoxinus dyads at ANP (r > 0.0, 1.5% of all possible dyads, Table 2), including two first-order relatives (r > 0.375) that were captured 2.17 km apart in different forest patches (the two patches were separated by ~1.5 km of matrix), suggesting a recent long-distance dispersal event across the grassland savannah matrix (Figure 2). All distantly related dyads were found further apart from each other (494 m - 2.5 km apart) than second-order relatives (up to 294 m). While the partners from 12% (2/17) of the second-order dyads were found in different forest patches, this was the case for 88% (7/8) of the distantly related dyads (Figure 2). No related (r > 0.0) E. myoxinus individuals were identified at MCF but sample size was rather limited (n = 14).

We identified 10 related dyads in R. rattus at ANP (r > 0.0, 1.1% of all possible dyads) and six at MCF (r > 0.0, 0.85% of all possible dyads, Table 2). All first and second-order relatives of R. rattus (r > 0.1875) were trapped within 110 m of each other, while distantly related R. rattus (r < 0.1875) were trapped within 760 m of each other. Related R. rattus were trapped in different forest patches only in one case (10%, 1/10, Figure 2).

In both E. myoxinus and R. rattus, second-order relatives were captured in closer proximity to each other than distantly related dyads (E. myoxinus: U = 16, p < 0.001 R. rattus: U = 28, p = 0.006). For both second-order and distant relatives, R. rattus dyads were in closer proximity to each other than the respective E. myoxinus dyads (second-order relatives: U = 136, p < 0.001, distant relatives: U = 26, p < 0.001).

We found significant effects of spatial autocorrelation with Mantel correlograms in E. myoxinus at ANP and for R. rattus at both ANP and MCF, but spatial autocorrelation could not be tested for E. myoxinus at MCF due to small sample size (Figure 3). For E. myoxinus at ANP, spatial autocorrelation between genetic distance and geographic distance was significantly positive for the first two distance classes (0 - 200 m, 200 - 400 m, Figure 3), while for R. rattus at ANP and at MCF, spatial autocorrelation was significantly positive only at the first distance class (0 - 200 m) (Figure 3).

3.3. Connectivity and Genetic Structure on the Landscape Level

Significant IBD effects were identified in both species (Figure 4). Genetic distance increased with increasing geographic distance in E. myoxinus at ANP (Mantel test, r = 0.222, p < 0.01, n= 52), and in R. rattus at both ANP (Mantel test, r = 0.4399, p < 0.01, n= 43) and MCF (Mantel test, r = 0.3211, p < 0.01, n= 37). No such effect was detected for E. myoxinus at MCF (Mantel test, r = -0.02, p = 0.52, n= 14) possibly due to small sample size.

Separate ADMIXTURE analyses for each landscape and species suggested an optimum of one or two clusters for each species within each landscape (see Figure S1 for justification and for results K = 3), of which we plotted the relative representation of K = 2 clusters in space (Figure 5). There was no strong signal of genetic structure in E. myoxinus in ANP, although one cluster (= blue) was dominant in the southwestern part of the landscape and rarer in the other parts (Figure 5a). At MCF, the southern and northern E. myoxinus subpopulations fell into largely separate clusters, with individuals showing only low proportions of admixture (Figure 5c). In contrast, there was no clear signal of genetic structure in R. rattus in either study region (Figure 5b,d).

At ANP, effective migration surfaces for E. myoxinus suggested highly variable connectivity across the landscape (Figure 6a). High effective migration rates (log(m) ≥ 0, = blue parts) were inferred within the central and northern sections of the fragmented landscape and in parts of the eastern continuous forest. In contrast, lower effective migration rates (log(m) < 0, = orange parts) were inferred for an area encircling the central part of high connectivity and for one stretch of the eastern continuous forest. At MCF, the few available samples came mainly from the two continuous forest sites in the North and in the South. We inferred high effective migration rates of E. myoxinus within each continuous forest site, while the savannah between the two forests acted as a barrier to gene flow (Figure 6c). At ANP for R. rattus, a reasonably high effective migration rate was inferred for large parts of the fragmented landscape and in the eastern continuous forest (Figure 6b). A low effective migration rate was only inferred between these two areas. At MCF, the effective migration surfaces showed a more complex picture for the northern and southern sites, respectively (Figure 6d). In the northern sub-section (6e) only a small number of individuals were present (n = 8) with mixed results; between two fragments there was a low effective migration rate, while between another set of fragments and the continuous forest there were high effective migration rates. In the southern forest sites (Figure 6f), there was a high effective migration rate amongst individuals in the continuous forest and along a section of the Mariarano river, but low rates were inferred amongst individuals in the forest fragments of the south and the west.

ResDisMapper also detected heterogeneous pattern of landscape resistance across the two study regions for the two species (Figure 7) with positive resistance values (green) indicating high levels of connectivity and negative resistance values (red) indicating low connectivity [63]. At ANP, areas of significantly high and low landscape resistance were inferred for E. myoxinus (Figure 7a). Overall, the Central and Southern fragments were all characterized by very low resistance, while very high resistance values were inferred towards and between the Northern fragments. At MCF, there were no significant resistance levels detected for E. myoxinus, but the individuals in the southern continuous forest had low resistance values while the individuals in the north had higher values (Figure 7c). At ANP for R. rattus, we detected significantly low resistance values in the western fragments and the eastern continuous forest, while significantly high resistance values were detected in the southern and central fragments (Figure 7b). At MCF, we detected significantly low resistance values for individuals in the southern continuous forest and along the Mariarano river, while we detected significantly high resistance values amongst the western fragments and in and among the northern forest patches (Figure 7d). In the northern sub-section (7e) there were only positive resistance values for the individuals within and between forest patches. In contrast, at the southern forest sites (Figure 7f), there were negative resistance values amongst individuals in the continuous forest and along one section of the Mariarano river, but also positive resistance values between fragmented forest sites and along other sections of the river.

3.4. Ecological Determinants of Genetic Diversity

Across all models and among all variables, region and species had the strongest impact on genetic diversity (Tables 3 and S5). Overall, E. myoxinus had significantly higher inbreeding coefficients (F_IS) and significantly lower observed and expected heterozygosity than R. rattus (Figure 8, Table 3). The inbreeding coefficient F_IS was significantly higher at MCF than in ANP in both species (Figure 8, Table 3 and Table 4). However, while observed and expected heterozygosity in E. myoxinus were higher at MCF than at ANP, the opposite regional effect was detected in R. rattus with higher values at ANP for both measures of heterozygosity (Figure 8, Table 4). Among all ecological variables, only the distance to edge had significant effects on genetic diversity. First, it was positively associated to the inbreeding coefficient in the overall sample (Table 3). Second, observed heterozygosity decreased significantly with increasing distance to edge in E. myoxinus, but not in R. rattus. No other parameter had any significant influence on genetic diversity in any model (Table 3 and Table 4).

4. Discussion

4.1. Do endemic and Native Rodents Use Non-Forest Matrices for Local Movement or Dispersal?

The spatial distribution of related individuals suggests that endemic E. myoxinus are able to move between habitat patches across matrix despite their primarily arboreal lifestyle [33]. Several first order and second order relatives were separated by savannah matrix at ANP, including one case of first order relatives who were found over two km apart from each other. The ability of E. myoxinus to cross matrix may explain why previous studies have found individuals in small fragments that have been isolated for at least 50 generations [21], although the ability to disperse to these fragments does not necessarily ensure long-term persistence in the fragments. The Mantel correlograms for E. myoxinus showed significant positive correlations at low distance classes followed by negative correlations at higher classes, suggesting some degree of local clustering of related individuals. However, both the steep drop-off of correlations after 400 m and the negative correlation values beyond suggests dispersal capacity is reduced in a fragmented landscape. When habitat is continuous, spatial correlation would be expected to decline in a more gradual pattern [69] and negative correlation values have been linked with reduced dispersal due to habitat disturbance [70]. Although the dispersal distances of E. myoxinus in continuous forest habitats are not yet known, our results suggest that despite reduced dispersal in fragmented landscapes, E. myoxinus maintains some connectivity between adjacent forest fragments and via opportunistic long-distance dispersal events across the matrix.

The invasive R. rattus showed different patterns of relatedness and dispersal across the landscapes. Related dyads were found in neighboring forest patches only in a few cases (2/16 cases), we did not observe any long-distance dispersal events, and all first and second-order relatives were captured within ~ 100 m of each other. These results suggest a high degree of philopatry and kin-biased spatial clustering in R. rattus which has already been described from other study sites in Madagascar [40]. This finding is also supported by the results of the Mantel correlograms as significantly positive spatial autocorrelation could only be shown at short distances (< 200 m). In a study on congener Rattus norvegicus in an urban landscape, a decline in correlations was also observed beyond 200 m but in a much more gradual pattern [57]. Correlograms with steep declines and significant negative values may again be interpreted as filtering effects of fragmentation as in the case of E. myoxinus [69,70]. Thus, dispersal of R. rattus in our study seemed to be reduced due to fragmentation as well, although dispersal across the matrix was possible but probably rare due to a preference for staying close to their natal site. While R. rattus are generally able to exploit a breadth of habitat types, our results indicate that they do not appear to use savannah as habitat, although it is clearly permeable to dispersal.

We found related R. rattus dyads in closer proximity than related E. myoxinus dyads. These differences can be attributed to differences in social organization and possibly resource requirements between the two species. E. myoxinus is a solitary, frugivorous, arboreal rodent [33,34] and was found in lower numbers in this fragmented landscape than R. rattus [21]. In contrast, R. rattus forms gregarious groups of related members and has an opportunistic omnivorous diet [40]. While the longer dispersal distances of E. myoxinus are in concordance with the solitary lifestyle of a dietary specialist [33], their lower reported abundance in fragments may rather be attributed to specific resource and microhabitat requirements than to a general lack of vagility between the continuous forest and/or forest fragments. In contrast, R. rattus seems to be largely philopatric with less dispersal across the matrix as a by-product of individuals preferring to stay near kin rather than an inability to disperse across open spaces. While our study did not directly test for this, it is also possible that human-mediated dispersal is occurring in these landscapes. Thus, while R. rattus and E. myoxinus both have different life history strategies, both species are able to maintain dispersal across non-habitat matrix.

4.2. Are Connectivity and Genetic Structure Affected Primarily by Isolation by Distance or by Habitat?

E. myoxinus at ANP showed a pattern of isolation by distance. Previous research on Eliurus spp. did not find significant effects of IBD using mitochondrial DNA [36,37], but we examined IBD at a smaller spatial scale (5 km vs 42 – 83 km) that is potentially affected by the small-scale dispersal patterns of individuals. The effect of scale can also be seen between the two landscapes, as the two E. myoxinus sub-populations at MCF showed strong signals of genetic differentiation between the two large forest patches separated by savannah, while a more complex patchwork of connectivity was observed at ANP not directly connected to a matrix-forest dichotomy. This effect is in line with range-wide patterns of E. myoxinus which has a wide distribution and is mainly structured by altitudinal gradients [32] with high rates of gene flow within the core of their range [34]. However, there has so far been little research done on this species at the scale of our study. The complex spatial patterns in connectivity at ANP could be due to other factors such as spatial variation in predation risk [71] or microhabitat preferences in E. myoxinus which may also explain the lack of unifying ecological patterns seen in other studies [16,21,72]

Rattus rattus sub-populations at both ANP and MCF showed significant effects of IBD. However, we also observed different patterns of connectivity across the landscape suggesting that other habitat factors such as matrix composition and quality may also structure these sub-populations. We observed two distinctive areas of high connectivity in larger forest patches at ANP and a more complicated pattern at MCF but similar areas of high connectivity in larger forest patches. These movement patterns suggest a general preference for larger patches of forest but also may be explained by the population dynamics of R. rattus in the landscapes. More specifically, the high degree of genetic connectivity among the larger fragments in ANP points towards a larger and well-established population nucleus consisting of multiple, largely philopatric matrilines that are rather isolated from neighboring lineages inhabiting the smaller, isolated fragments due to a scarcity of long-distance dispersal. Anthropogenic influences may also be structuring R. rattus populations in the region. R. rattus have been observed to follow human paths and settlements rather than dispersing in a Euclidean manner [42]. There are no long-term human settlements between our capture locations at ANP, although some villages are present north and northeast of the landscape [19], while the landscape at MCF is interspersed with several villages and agricultural land [45].

From a conservation management standpoint, our study showcases the difficulty in managing invasive R. rattus populations in Madagascar. Eradication is unfeasible, especially in a country with limited conservation resources. R. rattus are able to exploit and colonize continuous forest and forest fragments likewise, and the only habitat type they seem not to use, savannah matrix, cannot be expanded as it would also negatively affect endemic species such as E. myoxinus. It is likely that some of habitat-related factors structuring these two species are similar, such as vulnerability to predation [71] while microhabitat requirements may vary due to differing habitat preferences and feeding habits. While we found IBD is an important factor structuring these populations, these other habitat-related factors that are structuring populations should be examined further.

4.3. What Ecological Factors Affect the Genetic Diversity of Rodent Sub-Populations?

Only two ecological factors influenced the genetic diversity of both rodents and these were region, species and, in the case of E. myoxinus, distance to edge. The regional differences in the case of E. myoxinus were rather unexpected, since both the inbreeding coefficient as well as observed and expected heterozygosity were higher in MCF compared to ANP. These rather inconsistent findings may have been driven by the small sample size from MCF which precludes further interpretation. In congruence with the overall significant and positive relationship between F_IS and distance to edge, H_O decreased significantly with increasing distance to edge, but only in E. myoxinus. Taken together, these findings suggest that genetic diversity, at least in E. myoxinus, was higher closer to edges than towards the center of forests. Previous research from the same two landscapes found a decreased abundance of E. myoxinus in fragments with a higher percentage of edge zone, and it was concluded that E. myoxinus shows negative effects of habitat edges and of extensive habitat fragmentation [21]. If that is the case, a decrease of genetic diversity with increasing distance to edges may not be explained with a higher abundance of E. myoxinus near edges. Alternatively, the higher levels of genetic diversity near edges may have resulted from an increased diversity of mating opportunities with unrelated conspecifics near edges, since immigrants from other fragments arrive first in the edge zone of any given site, and long-distance dispersal was confirmed for E. myoxinus in the present study. Radespiel et al. [27] found decreased gene flow between mouse lemurs from edge and interior zones, suggesting a “local preference model” in which the ecological differences between edge and interior habitat coincide with a preference for the natal habitat. This may also be the case in E. myoxinus if individuals born in edge zones of fragments rather disperse across the matrix to other edge zones than towards interior sub-populations. Such a behaviour would lead to an increasing genetic diversity in edge zones despite a lack of high-quality habitat creating an ecological trap. Molecular edge effects are understudied and represent an important upcoming topic in fragmentation research, especially given the adoption of next-generation sequencing which may be able to parse out smaller-scale patterns in genetic differentiation [73]. Given that nearly half of the remaining forest in Madagascar is within 100 m of the edge [6], major conservation concerns would arise if forest edges are confirmed to act as an ecological trap for E. myoxinus or other endemic rodents.

The only significant factor affecting the genetic diversity of R. rattus when analyzed individually was region; we observed higher levels of genetic diversity at ANP. This may reflect the invasion dynamics of R. rattus in northwestern Madagascar. At least two founder events occurred in Madagascar, one in the North and one in the South, followed by a relatively slow range expansion across the island. The populations in northwestern Madagascar were suggested to have resulted from the southern invasion [42,74]. Consequently, the MCF population of R. rattus would lie at the edge of the range of the southern expansion and thus may have been more recently colonized than ANP and undergone a founder effect which is known to coincide with a population bottleneck and lower genetic diversity. The lack of other significant ecological impacts on genetic diversity may also suggest ecological flexibility, lack of habitat preference, and their principal ability to cross the matrix and mitigate its isolation effects. Future research, using methods such as approximate Bayesian computation (ABC) to infer demographic history of R. rattus could help tracing and reconstructing the original invasion(s) into Madagascar and its spread throughout the North into our study landscapes.

Genetic diversity of R. rattus was higher compared to E. myoxinus at both sites. Invasive species with a limited number of invasion points in recent history should have lower levels of genetic diversity compared with established populations of native species in a given landscape due to founder effects [75]. While previous research suggests R. rattus and E. myoxinus population abundances are not correlated in this site, suggesting a lack of competition [20,21], our results suggest that R. rattus are better equipped to expand into and establish viable populations within fragmented landscapes compared with E. myoxinus. Compared with R. rattus, E. myoxinus has more specialized habitat requirements, such as tree holes for sleeping sites and a frugivorous diet [35]. Thus, while direct competition is not necessarily occurring, the endemic E. myoxinus may still be at greater risk of extirpation in these fragmented landscapes over time, especially as deforestation continues and the amount of remaining habitat subject to edge effects increases. Given that distance to edges was the only significant ecological variable we found affecting E. myoxinus, we encourage further research in the microhabitat preference and natural history of this species in order to more effectively inform conservation.

5. Conclusions

We found variation in local movement, connectivity, and genetic diversity for rodents in Malagasy fragmented dry forests. Certain aspects of these species’ life histories seem to allow them to colonize and/or persist in fragmented landscapes. E. myoxinus, a solitary species with more specialized habitat requirements was identified to disperse over large distances between habitat patches across an open savannah matrix. However, decreased abundance, lower levels of genetic diversity overall, and lower genetic diversity with increasing distances from habitat edges suggest that the long-term persistence of this endemic species in fragmented landscapes is not secured. These overall threats are increased by the rapid rates of deforestation in Madagascar [6] and probably also affect other endangered local microendemic rodents across Madagascar for whom there is an even greater dearth of data. Conversely, R. rattus, a gregarious and ecologically flexible, invasive species, is well adapted to living in fragmented landscapes of northwestern Madagascar, although it maintains connectivity with short-distance rather than long-distance dispersal. Our results add further complexity to the literature on mammal populations in fragmented forests of Madagascar and suggest further research is needed. In particular, basic natural history information (e.g., dispersal patterns and distances, sleeping habits, social organisation, food resources, reproduction, microhabitat preferences) is missing for many species which would ultimately help to improve their effective conservation.