Evolutionary and demographic processes shaping geographic patterns of genetic diversity in a keystone species, the African forest elephant (Loxodonta cyclotis)

Abstract The past processes that have shaped geographic patterns of genetic diversity may be difficult to infer from current patterns. However, in species with sex differences in dispersal, differing phylogeographic patterns between mitochondrial (mt) and nuclear (nu) DNA may provide contrasting insights into past events. Forest elephants (Loxodonta cyclotis) were impacted by climate and habitat change during the Pleistocene, which likely shaped phylogeographic patterns in mitochondrial (mt) DNA that have persisted due to limited female dispersal. By contrast, the nuclear (nu) DNA phylogeography of forest elephants in Central Africa has not been determined. We therefore examined the population structure of Central African forest elephants by genotyping 94 individuals from six localities at 21 microsatellite loci. Between forest elephants in western and eastern Congolian forests, there was only modest genetic differentiation, a pattern highly discordant with that of mtDNA. Nuclear genetic patterns are consistent with isolation by distance. Alternatively, male‐mediated gene flow may have reduced the previous regional differentiation in Central Africa suggested by mtDNA patterns, which likely reflect forest fragmentation during the Pleistocene. In species like elephants, male‐mediated gene flow erases the nuclear genetic signatures of past climate and habitat changes, but these continue to persist as patterns in mtDNA because females do not disperse. Conservation implications of these results are discussed.

Field studies have shown strong evidence that, despite living in a fission-fusion society, female elephants remain with their closest kin after they mature . Genetic analyses have supported female philopatry by demonstrating almost complete uniformity of mtDNA haplotypes within families (Archie, Fitzpatrick, Moss, & Alberts, 2011). By contrast, male elephants upon reaching maturity disperse from their natal herds (Lee, Poole, Njiraini, Sayialel, & Moss, 2011) and enter periods of musth characterized by competitive interactions with other males for reproductive access to females . The dispersal of male elephants from their natal social groups thus mediates nuclear gene flow (Ishida et al., 2011;Nyakaana & Arctander, 1999;Roca et al., 2005Roca et al., , 2015. The phylogeographic patterns revealed by analyses of Y-chromosome sequences is similar to the pattern for other nuclear markers, but different from patterns shown by mtDNA, supporting the role of males in establishing nuclear phylogeographic patterns (Roca, Georgiadis, & O'Brien, 2007;Roca et al., 2005).
Because mitochondrial phylogeographic patterns are often discordant from nuclear patterns in species in which only males disperse (Petit & Excoffier, 2009;Toews & Brelsford, 2012), including elephants (Roca et al., 2005(Roca et al., , 2007, there is a strong need to analyze nuclear markers among forest elephants to examine more completely their evolutionary history and population structure. Furthermore, Central Africa, the region in which most forest elephants live, has suffered from the highest levels of elephant poaching of any subregion within the continent (CITES, 2012), and has been the main source for the illegal trade in elephant bushmeat and ivory . Forest elephant numbers declined by ca. 62% between 2002 and 2011, to <10% of their estimated historical population size, mainly due to illegal poaching for their tusks (Maisels et al., 2013). There is thus a strong need to examine fine-scale population substructure within forest elephants using nuclear markers, for proper conservation management of the species.
Here, we use microsatellite markers to examine nuclear genetic structure in the forest elephant. Ninety-three individuals from five localities in Central Africa and one individual from Sierra Leone were genotyped using 21 microsatellite markers. We examined nuclear genetic markers for geographic differences among forest elephant localities. We discuss the extent to which regional populations may or may not be genetically distinctive, and the implications of these findings for forest elephant conservation. We also specifically examine the degree of discordance between the phylogeographic patterns inferred using microsatellite markers and patterns previously reported for forest elephant mtDNA across the same tropical forest localities within Central Africa.

| Samples
This study was conducted under the University of Illinois Institutional Animal Care and Use Committee (IACUC)-approved protocol number 15053. Samples were collected in full compliance with required  (Olson et al., 2001) that historically included a mixture of forest and secondary grasslands suitable for both African elephant species (Groves & Grubb, 2000) (Tables 1 and S1).
As only one sample was available from SL, it was not included in some statistical analyses. Garamba has historically included mixed forest and savanna habitats suitable for both species of African elephant (Groves & Grubb, 2000). Most of our samples from Garamba are forest elephants, although a few are hybrids of savanna and forest elephants based on nuclear genotypes (Comstock et al., 2002;Groves & Grubb, 2000;Ishida et al., 2011;Roca et al., 2001;Rohland et al., 2010).  (Comstock et al., 2002;Groves & Grubb, 2000;Ishida et al., 2011;Roca et al., 2001;Rohland et al., 2010). Because our samples contain forest elephants from Garamba where a few hybrids of savanna and forest elephants have been identified based on nuclear genotypes (Comstock et al., 2002;Ishida et al., 2011;Mondol et al., 2015;Roca et al., 2001), it was determined that the microsatellites would be amplified in savanna elephants, as they would be needed to identify hybrids of savanna and forest elephants.
Details on the sampling and DNA extraction have been previously published (Ishida et al., 2011).

| Microsatellite genotyping
Allelic variation was examined at 21 microsatellite loci. These markers have been previously developed by Gugala et al. (Gugala, Ishida, Georgiadis, & Roca, 2016). Primer sequences are listed in Table S1. All forward primers included the M13 forward sequence (TGTAAAACGACGGCCAGT) at the 5′ end. The PCR primer mix consisted of a 5′ FAM-or VIC-fluorescent-labeled M13 forward primer (to label the PCR amplicon), along with the forward primer (with M13 forward sequence at the 5′ end) and reverse prim- Relying on a standard of known size, the binning function of the software Genemapper was used to determine fragment lengths, following the procedures indicated in the manual. For the DNA samples extracted from dung (BF), at least four independent amplifications were repeated to confirm homozygotes and three amplifications for heterozygotes (Allentoft et al., 2011).

| Characterization of microsatellites
Arlequin version 3.5.1.3 (Excoffier & Lischer, 2010) and GenAlEx 6.5 (Peakall & Smouse, 2012) were used to calculate expected heterozygosity (H e ) and observed heterozygosity (H o ). The software GenAlEx 6.5 was also used to calculate Shannon's diversity index (I) and to make allele frequency distribution histograms for each locus for each locality. Tests for Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) were conducted on the forest elephant microsatellite data. A Markov chain algorithm was used to test for HWE using 10,000 dememorization steps, 1,000 batches, and 10,000 iterations per batch using the software Genepop 4.2 (Rousset, 2008), and 1,000,000 steps and 100,000 dememorization steps were used with Arlequin version 3.5.1.3. LD was examined using 10,000 dememorization steps, 1,000 batches and 10,000 iterations per batch Pairwise F ST is shown for comparisons between localities using nuclear microsatellites (below diagonal) and mtDNA (above diagonal). F ST values of mtDNA are from Ishida et al. (2013). The values that are significant are indicated in boldface. Sample sizes are in parenthesis for mtDNA (top row) and microsatellites (first column). Localities corresponding to the abbreviations are shown in Figure 1. F ST values calculated using mtDNA are not shown for BF as comparable mtDNA sequences are not available for BF.

| Population genetic analyses
Analysis of molecular variation (AMOVA) was conducted in Arlequin version 3.5.1.3 (Excoffier & Lischer, 2010) using 10,000 permutations. For each pair of localities, Arlequin was also used to calculate pairwise F ST values and statistical support using 10,000 permutations, and to calculate the inbreeding coefficient for each locality, except SL, as only one sample was available from SL. STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000), which applies a model-based clustering algorithm to multilocus genotype data, was used to infer population structure using datasets that included or excluded the savanna elephants. STRUCTURE was run eight times for each value of K from 1-10, without the use of prior information on locality, under the admixture-correlated model, with each iteration using at least 1 million Markov chain Monte Carlo generations following a burn-in of at least 100,000 steps.
The uppermost hierarchical level of population structure was examined using an ad hoc statistic ∆K based on the rate of change in the log probability of the data for a given K between successive K values, implemented in Structure Harvester (Earl & vonHoldt, 2012).
Average coefficients were estimated for each K value that was estimated to be uppermost, and for lower values of K, employing the Greedy algorithm with 1,000 random input orders as implemented in the program CLUMPP version 1.1.2 (Jakobsson & Rosenberg, 2007). These outputs were visualized using DISTRUCT version 1.1 (Rosenberg, 2004). After identifying the hybrid elephants in GR, we also conducted STRUCTURE analyses excluding them to remove the influence of savanna elephant genotypes. We also conducted STRUCTURE analyses to compare each pair of localities.
Factorial correspondence analyses (FCA) were implemented in GENETIX version 4.05 (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 2004) to graphically plot the distribution of genetic variation for each locality with forest elephants (excluding the hybrids). Principal coordinate analyses (PCoAs) were implemented in GenAlEx 6.5 (Peakall & Smouse, 2012) to visualize the genetic relationships among individual elephants, both including and excluding the savanna elephants and hybrid GR elephants.

| Isolation by distance
The coordinate information of each locality was estimated using the LatLong.net website (http://www.latlong.net/), and the pairwise distances between each locality were calculated using the coordinate calculators and distance tools in GPS Visualizer (http://www. gpsvisualizer.com/calculators). To examine the relationship between genetic distances and geographic distances among forest elephants at the five localities in Central Africa, a spatial autocorrelation analysis was implemented in GenAlEx 6.5 (Peakall & Smouse, 2012).
The spatial autocorrelation analysis divided the pairwise distances into four ordinal classes and used 9,999 random permutations and 9,999 bootstrap iterations. Isolation-by-distance (IBD) analyses were conducted to test the relationships between genetic differences between each pair of localities in Central Africa and the geographic distance between them. These analyses used the Isolation By Distance Web Service Version 3.23 (Jensen, Bohonak, & Kelley, 2005). Two different measures of genetic distance were calculated: . Mantel tests were run with 30,000 randomizations (one-tailed test). For Slatkin's similarity index, we used the recommended log-transformation of both M and geographic distance. As we had only one sample from West Africa (from Sierra Leone, SL), we excluded this sample from analyses.

| Phylogenetic analyses
We inferred the phylogenetic relationships among localities using the neighbor-joining (NJ) method implemented in POPTREEW (Takezaki, Nei, & Tamura, 2014). As the number of samples was different among localities, we used D ST (Nei, 1972) and F ST (Latter, 1972) with sample size bias correction in addition to D A (Nei, Tajima, & Tateno, 1983) (for which sample size bias correction was not available) to calculate genetic distances among localities. Support for the nodes in each analysis was assessed using 10,000 bootstrap pseu-

| Characterization of forest elephants using microsatellites
Although 21 microsatellite markers were genotyped, three were removed before analyses. Marker Lcy-M4 had a low genotyping success rate due potentially to null alleles. In marker Lcy-M15, a 1-bp indel was detected in some savanna elephants; this marker was removed from the forest elephant analyses as some elephants in Garamba are savanna-forest elephant hybrids. The marker Lcy-M52 showed a significant deviation from HWE even after Bonferroni correction (p < .0026), and was monomorphic in three localities. The remaining 18 microsatellite loci (Table S2) Figure S1 and allele number, heterozygosity, and other information for each marker are listed in Table S2. The markers did not show high diversity in savanna elephants. This would be expected for two reasons. First, the markers had been designed based only on the presence of polymorphisms among forest elephants (Gugala et al., 2016) that are 4-7 million years divergent from savanna elephants (Brandt, Ishida, Georgiadis, & Roca, 2012;Rohland et al., 2010). Additionally, savanna elephants are known to have reduced nuclear genetic diversity relative to forest elephants (Roca et al., 2001;Rohland et al., 2010). In savanna elephants, allele numbers for the 18 microsatellite loci ranged from 1 to 4, with the mean of 2.39 ± 0.24. The mean H o was 0.21 ± 0.04, and the mean H e was 0.25 ± 0.05. Allele numbers, heterozygosity, and other information for each marker are listed for the savanna elephants in Table S3.

| Analyses of population structure across forest elephant localities
Population genetic analyses involved only forest elephants, except as otherwise indicated. Many analyses were run separately for datasets that included or excluded elephants in GR that were identified as hybrids between forest and savanna elephants, although the outcomes of these analyses were not greatly affected by including or excluding hybrids. Analysis of molecular variance (AMOVA) found that only 3.04% of the variance was accounted for by differences among localities (Table S4). F IT was 0.054 (p < .005) with significant deviation from HWE while F IS was 0.024 but did not deviate significantly from HWE. The value calculated for F ST was low at 0.030, and this value was statistically significant (p < .001).
F I G U R E 2 Bayesian clustering approach implemented in STRUCTURE (Pritchard et al., 2000) using 18 microsatellite genotypes, including both forest and savanna elephants or only forest elephants. (a) When both savanna and forest elephants were included, the uppermost K value was estimated as four (Earl & vonHoldt, 2012). At K = 2, the forest (Loxodonta cyclotis) and savanna (L. africana) elephants were almost completely separated into different partitions, with a few hybrids in GR, consistent with previous reports (Ishida et al., 2011). At higher levels of K, the additional partitions do not completely separate elephants from different localities, although forest elephants from the eastern Congolian forest (BF, GR) show different patterns in their partitions than localities further west. (b) When only forest elephants were analyzed, the uppermost K value was estimated as three. Partitioning within the forest elephants resembles that seen in panel A, with a distinctive but incomplete pattern of partitioning between eastern and western localities. Note that additional analyses ( Figure S2) excluding the forestsavanna hybrid elephants from GR did not greatly affect the patterns of partitioning

Loxodonta cyclotis L. africana
For each pair of forest localities, F ST was also calculated (Tables 1 and S5) for pairs of localities (excluding SL for which the sample size was one). We identified five elephants in Garamba as hybrids using STRUCTURE (see below) and removed these five elephants from the analyses. Pairwise F ST values were high when each forest locality was compared to savanna elephants, which for these analyses were grouped together (Laf in Table S5). In the comparisons involving a forest locality and the grouped savanna elephants, pairwise F ST ranged from 0.40 to 0.65, with a statistically significant result for each comparison. F ST was also calculated between each pair of localities containing forest elephants (Tables 1 and S5), with values ranging from zero (BF and GR) to 0.07 (GR and OD). F ST values were low between pairs of forest elephant localities although there were modest but statistically significant differences between some localities in the eastern and western regions of the Congolian forest block (Tables 1 and S5 Figure S6). Although the distant placement of the Sierra Leone elephant is intriguing, we caution that no strong conclusion can be drawn from a single individual  Figure   S2). We also conducted STRUCTURE analyses for each pair of localities. We detected differences but incomplete partitioning between GR and the two western localities, DS and LO in the pairwise comparisons ( Figure S3). Different patterns of partitioning were not observed between BF and other localities in these pairwise comparisons, presumably due to the small sample size for BF ( Figure S3).  (Figures 3b, S4b). This tendency was not evident using biascorrected D ST and F ST (Figures 3b, S4b). In all three trees, savanna elephants (Laf) formed a lineage distinct from forest elephants ( Figure S4b). Among forest elephants, the eastern Congolian forest localities (DS, OD, and LO) formed a clade that was distinct from a clade formed by the western Congolian forest localities (BF and GR), while SL (in the Guinean forest block) had a long branch separating it from all other forest elephant localities (Figures 3b, S4b). To examine whether the separation of SL on the tree was merely due to it having the smallest sample size of n = 1, we reran the analyses while also limiting the sample size of LO to a single individual ( Figure   S6). Reducing the sample size of LO to one individual affected the trees based on D A , with the terminal branch length of LO becoming relatively longer. By contrast, trees based on D ST and F ST that were corrected for sample size bias consistently showed a relatively long branch for SL but not for LO when a single individual was used for each locality ( Figure S6). This was true for three different individuals from LO, randomly chosen and used in separate analyses in which the sample size of LO was limited to one. For the two bias-corrected methods, the long branch length was consistently present for SL with n = 1, but not for LO with n = 1, which may suggest that the long separation between SL and the other populations may be due to actual genetic differences, and not be a mere artifact of the small sample size.

| Evidence for isolation by distance among forest elephants in the Congolian forest block
We examined the degree to which genetic differences among forest elephants at different localities varied with the geographic distances separating them. Geographic distances were computed between each pair of elephants, with the distances placed into quartiles (x-axis in Figure 4a) to implement a spatial autocorrelation analysis (Peakall & Smouse, 2012). Genetic distances between pairs of elephants were also determined (y-axis in Figure 4a). A spatial autocorrelation analysis showed a correlation between genetic distance and geographic distance. Forest elephants that were close to each other geographically were also more similar genetically than were the geographically distant forest elephants ( Figure 4a). Support for isolation by distance (Jensen et al., 2005) (Figure 4c).

| Current population structure among forest elephants
Nuclear genetic differences were detected among Central African Hybrids between forest and savanna elephants have been documented, but only within relatively narrow transition zones between forest and savanna habitats (Comstock et al., 2002;Ishida et al., 2011;Mondol et al., 2015;Roca et al., 2001). here with those previously developed, is likely to provide greater F I G U R E 4 The relationship between genetic distance and geographic distance among forest elephants from different localities. (a) Spatial autocorrelation analysis found that genetic distance was correlated with geographic distance (r is the spatial autocorrelation coefficient; U is the upper 95% randomization limits of r; L is the lower 95% randomization limits of r). (b) Results consistent with potential isolation by distance were obtained by comparing genetic distance (F ST ) to geographic distance (km) in pairwise comparisons of forest elephant localities (r = .86, p < .01). (c) Results consistent with potential isolation by distance were obtained when genetic distances were calculated using Rousset's distance F ST /(1-F ST ) and compared to geographic distances (km) in pairwise comparisons of the genotypes of forest elephants grouped by locality (r = .85, p < .01). Five forest-savanna elephant hybrids from GR were excluded from the analysis shown in each panel; Sierra Leone (SL) was not used in the pairwise comparisons because the sample size was one. For results including the SL sample, see Figure S7  Geographic distance precision in estimating the degree to which hybrids received alleles from one lineage or the other (Boecklen & Howard, 1997). It would also be likely to further increase the accuracy and precision of estimating the provenance of confiscated ivory using nuclear markers .
The distinctiveness of elephants from West Africa has been proposed (Eggert et al., 2002), but based largely on mitochondrial DNA data, which can be misleading in elephants (due to maternal inheritance and female philopatry) (Ishida et al., 2011(Ishida et al., , 2013. The single individual from West Africa (Sierra Leone) in this study appeared to anchor one of the axes in the FCA analysis, and also generated a long branch in the phylogeny, when compared to other isolated individuals (Figures 3 and S4). While suggestive, this is not conclusive evidence for the distinctiveness of West African elephants (a larger sample set is required). However, the Benin/Dahomey Gap separating the Congolian from Guinean forest blocks may hinder gene flow in forest elephants, as it has affected the distribution of other forest-dwelling taxa (Linder, 2014). For this reason, we previously recommended that a conservative approach would treat elephants on either side of the Gap as deserving of separate conservation status (Roca et al., 2015). Whether Guinean and Congolian forest elephants form genetically distinctive groups (based on nuclear DNA analysis) remains one of the most important unanswered questions in elephant conservation genetics (Roca et al., 2015).

| Discordant mitonuclear patterns and the role of range expansion
Mitochondrial DNA patterns have previously been examined in African forest elephants (Debruyne, 2005;Debruyne et al., 2003;Eggert et al., 2002;Ishida et al., 2013;Johnson et al., 2007;Nyakaana et al., 2002), revealing five mitochondrial subclades with distinctive geographically restricted distributions (Ishida et al., 2013). In a previous study (Ishida et al., 2013), pairwise F ST values calculated using mtDNA were quite high, ranging from 0.38 to 0.87 when estimated pairwise between localities (Table 1). However, in taxa with malebiased dispersal, phylogeographic patterns are often discordant between nuclear and mitochondrial DNA (Petit & Excoffier, 2009;Toews & Brelsford, 2012), and discordant mitonuclear patterns have been reported among living and extinct elephantid species (Enk et al., 2011;Lei, Brenneman, Schmitt, & Louis, 2012;Meyer et al., 2017;Palkopoulou et al., 2015;Roca, 2015;Roca et al., 2005). This is consistent with the current analysis, which found nuclear genetic differentiation among forest elephants to be much lower than what might be inferred using mtDNA alone, with all values of F ST ≤ 0.07 among forest elephant populations across Central Africa (Table 1).
We would note that the high mutation rate of mtDNA would not account for the discrepant patterns, because a faster-evolving marker would only reveal greater resolution than a slower one; by contrast mtDNA demonstrates a strikingly different phylogeographic pattern than nuclear DNA in forest elephants (Figure 2 vs. Figure S8) as in other elephantids (Enk et al., 2011;Lei et al., 2012;Meyer et al., 2017;Palkopoulou et al., 2015;Roca et al., 2005Roca et al., , 2007Roca et al., , 2015.
Mitonuclear discordant patterns in most cases have been attributed to adaptive introgression of mtDNA, demographic disparities and sex-biased asymmetries, with some studies also implicating habitat changes and hybrid zone movements (Toews & Brelsford, 2012). In the case of forest elephants, adaptive introgression of mtDNA is unlikely, not only because selective sweeps are unlikely to occur in markers carried only by the nondispersing sex (Petit & Excoffier, 2009), but because mtDNA shows greater differentiation across the forest elephant range than does nuDNA (Table 1). Instead, sex-based differences in gene flow appear to be responsible for the discordant mitonuclear patterns, with some impact likely due to changes in habitat across geological time. Because female elephants are matrilocal and remain with their natal social group Hollister-Smith et al., 2007), this behavior can account for the persistence of geographic structuring in forest elephant mtDNA (Ishida et al., 2013). By contrast, male elephants disperse from their natal social groups and mediate nuclear gene flow across the landscape (Nyakaana & Arctander, 1999;Roca et al., 2005Roca et al., , 2015. In forest elephants, the mitonuclear patterns were likely impacted by habitat changes across geological time, and discordant mitonuclear patterns may provide a means for studying their range expansion after the end of the last glacial period. Genetic patterns largely depend on the demographic and ecological characteristics of a species (Castric & Bernatchez, 2003). Spatial patterns of genetic diversity may also reflect past changes in climate and habitats that expanded and contracted the ranges of species, sometimes at a fast pace (Hewitt, 2000). Current spatial genetic diversity may reflect such past events rather than species demography, with geographic differences in genetic diversity reflecting the effects of past climate-driven range dynamics (Hewitt, 2000). Range expansion can lead to patchiness after migration, due to long-distance movements followed by population expansions (i.e., a leptokurtic distribution of dispersal distances during colonization) (Ibrahim, Nichols, & Hewitt, 1996;Klopfstein, Currat, & Excoffier, 2006;McInerny, Turner, Wong, Travis, & Benton, 2009). These compounded foundation processes can lead to increased genetic differentiation, although such effects are negatively correlated with migration rate, because migration decreases lags in colonization and reduces the strength of founder effects (Klopfstein et al., 2006;McInerny et al., 2009).
In many species, it is often difficult or impossible to infer processes that are not directly observable from the current spatial genetic structure, especially as various processes may create similar patterns (McIntire & Fajardo, 2009). However, in elephants extreme sex differences in dispersal may allow for the study of both current demographic effects through the examination of male-mediated nuclear patterns (this study), and for the study of ancient landscape effects through analysis of mitochondrial patterns mediated by females ( Figure S8) (Ishida et al., 2013), which may retain signatures of leptokurtic dispersal and compounded foundation processes.
Pleistocene glacial cycles caused habitat changes that temporarily isolated populations of some species, with repeated cycles of isolation followed by expansion and contraction (Hewitt, 2000). During periods of spatial expansion, alleles present at the expanding edge of the species range can reach high frequencies (Klopfstein et al., 2006), with expanding populations potentially subject to iterated founder effects (Klopfstein et al., 2006). During expansions, rare long-distance dispersal events followed by exponential population growth can generate long-term patchiness in population structure (Hewitt, 2000;Ibrahim et al., 1996). Such ancient events may be preserved in elephant mitochondrial geographic patterns, which are likely to be stable due to female philopatry.
Further, studies of forest elephants would also avoid a common pitfall of using too small a geographic scale (Jenkins et al., 2010) while benefiting from a very large number of museum samples that are available for mtDNA analyses and also have precise provenance information. Species refugia during glacial cycles are better characterized in Europe and North America than elsewhere, and the forest elephant may provide novel insights into the impact of global glacial cycles in the African tropics (Hewitt, 2000).

| The potential role of isolation by distance
Distinguishing between discrete population genetic structure (Figures 2 and 3) and isolation by distance ( Figure 4) can be difficult (Meirmans, 2012). Models of isolation by distance are often used to approach the balance between drift and dispersal (Jenkins et al., 2010;Wright, 1943). Limited migration permits genetic drift, increasing population differentiation and leading to a correlation between neutral genetic differentiation and geographic distance (Jenkins et al., 2010;Wright, 1943). IBD develops in continuously distributed species when divergence accumulates due to genetic drift between locations separated by geographic distances large enough to overcome the homogenizing effects of gene flow (Chesser, 1983;Crispo & Hendry, 2005;Hewitt, 2000;Jenkins et al., 2010;Wright, 1943), and is common in natural populations (Jenkins et al., 2010). IBD can remain at equilibrium, after sufficient time has elapsed for genetic patterns to be established and stabilized (Castric & Bernatchez, 2003).

Forest elephants have been contiguously distributed across
Central Africa from the start of the Holocene (Plana, 2004), until at least the mid-to late 1900s (Douglas-Hamilton, 1987). The correlation between genetic and geographic distances among localities in Central Africa (Figure 4) may be attributable to isolation by distance (IBD) (Jenkins et al., 2010;Wright, 1943). However, because forest elephant populations expanded from discrete glacial refugia, it would be difficult to distinguish a role of IBD from the persistence of discrete population structure (Meirmans, 2012) that would have been diminished but perhaps not erased by postglacial gene flow.

| Conservation implications
Given the endangered status of forest elephants, and their role as a keystone species, discussion of the conservation implications of our results is warranted. An important step for the conservation of forest elephants would be universal recognition of its status as a separate species in need of species-specific conservation measures (Roca et al., 2015). Central Africa, the region in which most forest elephants live, has suffered from the highest levels of elephant poaching of any subregion (CITES, 2012) and has been the main source for the illegal trade in elephant bushmeat and ivory  leading to massive declines in their numbers (Maisels et al., 2013).
The increase in human numbers and activities has caused fragmentation of elephant habitats and range (Blake et al., 2007(Blake et al., , 2008. Such fragmentation reduces genetic connectivity, which can lead to loss of genetic variation, and ultimately to inbreeding and increased drift. Recognition both of species divisions and of population genetic patterns below the species level is essential for the maintenance of biodiversity, and an important conservation principle is to retain populations representing existing genetic variation (Moritz, 1994;Ryder, 1986). Our findings suggest that West African, western Congolian, and eastern Congolian forest elephants should be managed separately. For species such as the forest elephant that exhibit patterns indicative of limited population structure (Figures 2 and   3) and possibly isolation by distance (Figure 4), Chesser (Chesser, 1983) has suggested dividing the range into management units, with greater genetic exchanges within units than across units in order to balance the need for connectivity with the need to prevent loss of alleles. Within regions, it is important to prevent extreme habitat fragmentation and retain connectivity so that gene flow can continue among populations (Chesser, 1983). Should the destruction of elephants cease while habitats remain, populations could expand from multiple locations. Should active management become necessary for recolonization, the best source populations for translocations would be those that are geographically close, as these would be most similar to any extirpated population (Monsen & Blouin, 2004).
Forest elephants play a critical role in shaping their ecosystem, maintaining tree diversity, dispersing seeds in greater quantities and distances than most other fauna, and improving rates of seed germination following passage through the gut (Campos-Arceiz & Blake, 2011). In the central Congo Basin, seventy-eight percent of the larger tree species in the rain forest are dependent on forest elephants for seed dispersal (Blake et al., 2007). Thus, it is important to consider their ecological role in seed germination and dispersal when determining conservation priorities. Terrestrial ecoregions of the world have been mapped to identify areas of high biodiversity and representative communities (Olson et al., 2001). The range encompassed by Central African forest elephants includes five tropical forest ecoregions, three ecoregions of forest-savanna mosaic, as well as the Albertine Rift montane forests and Cameroonian highlands forests ( Figure S9) (Olson et al., 2001). Given limited genetic differentiation among elephants across the Congo Basin, these ecoregions could serve as management units for the forest elephants, reflecting the dependence of many plant species on elephants for seed germination and dispersal. Within regions, the extirpation of local elephant populations should be minimized, to minimize effects on the regional flora.
An important further consideration is the impact of forest elephant conservation on those plant species that are rare or regionally restricted. A survey of 5,881 species of plants across sub-Saharan Africa has been used to map endemism richness, which is defined as the sum of species present at a geographic location, but with the occurrence of each plant species inversely weighted by the size of its range ( Figure S10) (Linder, 2014;Linder et al., 2012). The endemic richness of plants has been found to be highly congruent with the endemic richness of frogs, snakes, birds, and mammals, which when combined indicated that the Benin Gap has influenced these patterns (Linder, 2014).
Congruence across these groups has been attributed to (1) a similar influence across vertebrates of the vegetation and flora, (2) common responses to the same climatic parameters, and (3) a common underlying history (Linder, 2014). For plants, regions with highest levels of endemism have been identified (Linder, 2014;Linder et al., 2012) and are indicated in Figure S10. These regions of endemism may be considered when setting conservation priorities for forest elephants, although a more precise mapping of endemism among those plants that are dependent on elephants would be helpful. Conserving elephants in the regions rich in plant endemism would directly benefit the conservation of many plant species dependent on elephants, and indirectly benefit other vertebrate and nondependent plant species for which levels of endemism are geographically congruent.

The work was funded by USFWS African Elephant Conservation
Fund Grants AFE-0778-F12AP01143 and AFE1606-F16AP00909.
The study was conducted under the University of Illinois Institutional Animal Care and Use Committee (IACUC) approved protocol number 12040. Samples were imported using appropriate CITES permits.
For technical or other assistance, we thank M. Malasky, R. Hanson,