Abstract

To accurately estimate the genetic diversity and population structure for improved conservation planning of Milicia excelsa tree, 212 individuals from twelve population samples covering the species' range in Benin were surveyed at seven specific microsatellite DNA loci. All loci were variable, with the mean number of alleles per locus ranging from 5.86 to 7.69. Considerable genetic variability was detected for all populations at the seven loci ( ; ). Moderate but statistically significant genetic differentiation was found among populations considering both (0.112) and (0.342). All of the populations showed heterozygosity deficits in test of Hardy-Weinberg Equilibrium and significantly positive values due to inbreeding occurring in the species. Pairwise values were positively and significantly correlated with geographical distances ( ; , Mantel's test) indicating that populations are differentiated by “isolation by distance.” Bayesian analysis of population structure showed division of the genetic variation into four clusters revealing the existence of heterogeneity in population genetic structure. Altogether, these results indicate that genetic variation in Milicia excelsa is geographically structured. Information gained from this study also emphasized the need for in situ conservation of the relict populations and establishment of gene flow corridors through agroforestry systems for interconnecting these remnant populations.

1. Introduction

Tropical trees suffered from several threats which have considerable long-term effects on the demographic and genetic structure of populations worldwide [1]. Deforestation, habitat fragmentation, and selective logging result in the loss of genetic diversity, the extinction of local populations, the reduction of population size, and disruption of mutualisms within pollinators and seed-dispersal agents [14] Lemes et al., 2003 [5]. To evaluate and reduce the genetic effect of deforestation and logging, it has become a priority to obtain information on the natural levels and distribution of genetic variation in population of tropical trees. Investigations on population genetic diversity and structure of populations within a species may not only illustrate evolutionary histories, processes, and mechanisms but also provide useful information for the biological conservation of endangered species [6].

Earliest studies of tropical trees have used isozymes as the primary genetic markers and have shown that most of the species investigated are outcrossed, exhibit high levels of genetic diversity and gene flow, and carry much of the variation within rather than among populations [2, 7]. More accurately, recent studies using microsatellites with strong discriminatory power have reported high level of inbreeding, restricted gene flow and isolation by distance in some tropical trees due to fragmentation [1, 8, 9] Lemes et al. [5].

From the early 1980s up to late 1990s, Milicia excelsa, (Moraceae) a large deciduous, dioecious species up to 30–50 m height, with a diameter of 1.70–2 m, commonly called iroko, was the focus of an important trade between Africa and Europe. Indeed, many African countries were exporters of iroko timber, especially Ghana (traded together with M. regia), Côte d'Ivoire, Cameroon, Congo, and Gabon [10]. In Benin, as the species ranks the first among valuable timber species, superior individuals were almost all logged from forests. Thus as a result of its overexploitation, the species was classified on the IUCN red list as close to vulnerable [11] and becomes the focus of conservation concern in many African countries. Therefore, there is an urgent need for effective conservation and management of the remnant populations. For this purpose, estimates of population genetic parameters are essential. The variability observed at different loci with codominant markers provides estimates of inbreeding, heterozygosity, gene flow, and outcrossing rate, all of which are important measures for assessing the conservation and management status of tropical trees under intense human pressure (Lemes et al. [5]).

In this work, we assessed the genetic diversity and population structure of the remnant natural populations of Milicia excelsa across its range in Benin using microsatellite markers recently developed for the species [12]. The specific questions addressed in this study were (i) what is the amount of genetic diversity and differentiation harboured within-and among populations? (ii) what are the implications for enhanced conservation of this species?

2. Materials and Methods

2.1. Plant Material

Milicia excelsa is a dioecious species. Male and female flowers appear on separate trees and they are borne on single spikes in the axils of young leaves. Male flowers are white in a slender catkin with 6 to 8 mm large, closely crowded on pendulous, slender catkins (spikes) 15–20 cm long, dangling from twigs of the outer crown. Female trees produce erect flower spikes about 5-6 cm long and 2 cm thick. Female flowers are greenish, in shorter and much fatter spikes, the styles of each flower project so that the inflorescence appears hairy. Fruits, 5 to 7.5 cm long with a diameter of 2 to 2.5 cm, are arranged along a longitudinal axis with 1 seed on each side, green, wrinkled, fleshy, and resembling a fat green caterpillar; no change in the colour of the syncarp when mature, but the flesh between the actual fruit softens. Seeds are hard, small, and lie in the pulp.

The fruit bat Eidolon helvum feeds on Milicia fruits and is considered as the principal agent of dispersal of this giant African tree both quantitatively and qualitatively. A part from bat, it was noticed that several other agents are involved in the Milicia seeds dispersion among which, the red-crowned parrot which defecate seeds primarily in the seed-hostile environment below parent tree crowns, elephants and monkeys or duikers which probably kill and defecate many seeds in large clumps which predisposes them to seed predation, seedling competition, or fungal attack. It is noteworthy that these agents are poor dispersers compare to fruit bat.

Milicia excelsa species is found in transitional vegetation between closed forest and savanna. It is often found in deciduous forests, semideciduous forests, and evergreen forest and sometimes in savannah woodlands. Occasionally it is found in isolated relict forest up to about 1300 m from sea level.

Young leaves samples from a total of 212 Milicia excelsa mature trees (15 to 21 individuals per population) were collected in twelve natural populations representing the geographical distribution of the species in Benin (see Table 1 for more details). Populations were natural, geographically distant for more than 50 km and trees were sparsely distributed within them. Trees sampled within a population were distant from at least 100 m. Leaves were silicagel-dried and kept in freezer at C until DNA extraction.

2.2. DNA Extraction and PCR Protocols

Total genomic DNA was extracted from about 0.05 g of dried leaves using the DNeasy Plant Mini Kit (QIAgen, UK, Ltd) and following the protocol provided by the manufacturer. DNA concentration was determined by electrophoresis on agarose-TAE gel stained with ethidium bromide [13].

Molecular variation at 7 microsatellite loci was investigated using the following primers designed for Milicia excelsa [12]: Mex163, Mex63, Mex137, Mex51, Mex202, Mex81, Mex95. Polymerase chain reaction (PCR) amplifications were carried out using a 96-wells thermal cycler PTC Peltier (MJ Research, Waltham, Massachusetts) and a 10- L reaction volume, with 10 ng DNA, 1 PCR buffer (Invitrogen, Carlsbad, California), 1.5 mM MgCl2 (Invitrogen), 0.4 nmol dNTPs (Invitrogen), 2 nM each of forward and reverse primers and 0.6 U of Platinum Taq polymerase (Invitrogen). 0.2 nM of fluorescent M13 ( CACGACGTTGTAAAACGAC ) labeled with IRDye700, was added to the end of forward primers. This technique is time consume less and cheaper than the use of fluorescent dCTP The PCR conditions were: initial denaturation at C for 4 min followed by 30 cycles of denaturation at C for 30 s, annealing at C for 30 s and extension at C for 3 min, with final extension at C for 15 min. [12]. PCR products were analyzed on 6.5% denaturing polyacrylamide gels (Long Ranger, Biowhittaker Molecular Applications, Rockland, Maine) using an IR2DNA analyzer (Li-cor). Allele sizes were determined using the v.2.0 software (Li-cor).

3. Numerical Analysis

3.1. Microsatellite Diversity and Population Differentiation

Test for deviation from Hardy-weinberg equilibrium (HWE) at each locus and genotypic linkage disequilibrium among seven pairs of microsatellite loci in each population were performed using Fisher’s exact tests with unbiased -values available in the GDA program version1.1 [15]. -values were corrected for multiple comparisons by applying a sequential Bonferroni correction [16]. Using the GDA program, genetic parameters for each population and overall loci were assessed by computing the average number of alleles per locus (AR) adjusted for the differences in sample size by resampling [17], observed heterozygosity ( ), and heterozygosity expected under Hardy-Weinberg equilibrium ( ). Nei’s estimates of heterozygosity for all populations at each locus including the observed heterozygosity ( ), the total genetic diversity pooled over all populations ( ), the mean genetic diversity within populations ( ) and the proportion of the total genetic diversity that occurs among populations due to different genotypes were calculated using FSTAT version 2.9.3.2 [18], available online at http://www2.unil.ch/popgen/softwares/fstat.htm. For purpose of comparison, genetic differentiation was quantified using both statistics and differentiation based on allele size ( ). Weir and Cockerham [19]; Excoffier et al. [20]; and Weir [21] estimators of , , , and were estimated for each allele, locus and overall using FSTAT version 2.9.3.2 [18]. Jackknifing over samples and loci and bootstrapping over loci were automatically performed for the statistics , , . Values of overall loci were estimated as weighted and unweighted [22, 23].

3.2. Geographic Structure

Geographical trends in the distribution of genetic diversity were investigated by three different methods. First, the null hypothesis of independence between pairwise and geographical distance was tested using Mantel test [24] implemented in ARLEQUIN program version 3.1 [25]. In addition, the Bayesian methods for the analysis of population genetic structure were used to estimate hidden population substructure by clustering sampled populations into panmictic groups by the Bayesian Analysis of Population Structure (program BAPS, [26], available online at http://www.rni.helsinki.fi/~mjs/). This does not require conditioning on known population structure. BAPS was run five times for 105 iterations after a burn-in period of 20 000, randomly mixing the order of populations in the input file. The resulting partitions were presented on a UPGMA dendrogram using Kullback-Leibler divergences among populations.

The model-based clustering algorithm was applied to identify clusters of genetically similar individuals and to test the proportion of genetic admixture among the clusters using STRUCTURE version 2.1 [27], available online at http://pritch.bsd.uchicago.edu/software/structure2_1.html. This was to assign individual multilocus genotypes probabilistically to a user defined number K of cluster or gene pools [28], achieving linkage equilibrium within cluster. We used an admixture model in which the fraction of ancestry from each cluster is estimated for each individual. A burn-in period of 1 000 000 iterations and 100 000 Markov Chain Monte Carlo (MCMC) repetitions were used and the program was run three times for each on the total dataset without any prior information on the origin of each sampled individuals and populations. Parameter of individual admixture alpha was chosen to be the same for all clusters and was given a uniform prior. The allele frequencies were kept independent among clusters to avoid overestimating the number of clusters [28, 29]. Proportions of ancestry were averaged over individuals within each population sample and the corresponding pie charts were plotted onto the Benin map. The clusters are referred to as gene pools [28].

4. Results

4.1. Within and among Population Variation at the Microsatellite Level

Almost all loci met Hardy-Weinberg expectations except Mex81 in Bante population and Mex95 in Tamarou population which showed significant deviation from HWE equilibrium after Bonferroni corrections . Exact test for genotypic linkage disequilibrium showed significant deviations for 7 out of 252 comparisons , which may be due to chance alone and these deviations are all associated with Mex81 and Mex95 suggesting that this disequilibrium could largely be due to the Hardy-Weinberg disequilibrium in these two loci since pairwise measures include both within- and among-loci disequilibrium. A total number of 147 alleles were recorded over the 7 loci in the 212 individuals with a range from 13 alleles (Mex202) to 32 alleles (Mex63) and an overall average of 6.805 Table 2. Most of the loci were variables across populations with mean observed heterozygosity ranging from Mex95 to Mex163 . Total gene diversity is high and similar for each locus with a range from Mex202 to Mex63 ; and gene diversity within populations ranging from Mex202 to Mex81 .

Allelic richness in populations varies from 3.71 (Sakete) to 5.86 (Lokossa) with an average number of alleles per locus within populations ( ) of 4.60 (Table 1). The number of private alleles (Pa) across populations ranged from 1 to 5 with a total number of 29. The mean percentage of polymorphic loci within populations ( ) was 100% (0.99 criterions). At the population level, the observed heterozygosity ranged from (Niaouli and Aplahoue populations) to (Save population) with an average value of 0.54. The average gene diversity within population was (a range from 0.732 to 0.855). Populations were on average inbred with a mean inbreeding coefficient of , indicating a significant excess of homozygotes relative to Hardy-Weinberg expectations (Tables 1 and 2). Bante and Tamarou Populations showed significant pairwise disequilibrium between Mex81 and Mex51, and Mex95 and Mex51 respectively .

4.2. Geographical Patterns of Population Genetic Structure

Genetic differentiation among populations estimated over loci by was moderate but within the range of the estimated values using (0.112 0.018). Mean values of and were 0.342 and 0.112 respectively suggesting that about 34% and 11% of the total variation resided among populations under stepwise-mutation model (SMM) and infinite allele model (IAM). Differentiation among populations based on the allele sizes, was significantly larger than differentiation based on allele identities , indicating that stepwise-like mutations contributed to overall among population differentiation. Based on the private allele method Gene flow estimated for the studied populations from the formula [14] equaled to migrant per generation suggesting a restricted extent of gene flow among populations.

Correlation between geographical and genetic distances drawn from Mantel test was significantly positive ( ; ) as shown by Figure 1, indicating an evidence for isolation by distance, that geographically close populations tended to be genetically more similar.

BAPS analysis showed division of the populations into four clusters almost identical to the Neighbor-Joining dendrograms based on Nei’s unbiased genetic distance. This result clearly showed the existence of heterogeneity in population genetic structure of Milicia excelsa. Southern east populations clustered to form group 1, southern west populations composed group 2, central and northern east populations grouped together except Biro which clustered with the first group. Population of Djougou remained alone in the fourth group. The model-based clustering method with STRUCTURE confirmed the heterogeneity in patterns of population genetic structure in Milicia excelsa in Benin although the inference of the number of gene pool K was not straightforward. Seeing that the log-likelihood values, ln Pr X/K increased with K values increasing, we choose the smallest value that captures the major structure in the data as suggested by Pritchard and Wen [27]. The value of brought us close to the results from BAPS analysis moreover ln PrX/K values vary among different runs for the same K value when (Figure 2). The proportion of an individual genome from each population that contributed to each of the four clusters under a model with the highest probability is given by Figure 3. Gene pool 1 is most abundant across central region’s populations; gene pool 2 is predominant in southern east populations while gene pool 4 mainly occurs in southern west populations. Gene pool 3 makes the bulk part of northern west population but is also remarkably represented in some populations from the south. Biro population which, clusters with southern east populations, is geographically located in the northern east and genetically shares some common features (gene pool 1, gene pool 2, and gene pool 3) with central, southern east and north populations, respectively.

5. Discussion

5.1. Microsatellite Diversity in Milicia excelsa

The overall pattern of genetic variation at microsatellite loci in Milicia excelsa revealed considerable genetic diversity ( , , and ). The very high level of polymorphism detected could be due to high number of repeats (14 to 29 repeats) in the used loci [12] as SSRs with many repeats have been shown to be more polymorphic than ones with few repeats [4, 31, 32]. Genetic diversities observed in Milicia excelsa (AR, HO, HE and HS) are lower than the pattern found in other microsatellite studies of tropical tree species including Symphonia globulifera [1], Carapa guianensis [4], Swietenia humilis [8], Caryacar brasilense [9], and Swietenia macrophylla (Lemes et al. [5]) but higher than genetic diversity observed in Vitellaria paradoxa [33].

5.1.1. Population Structure and Inbreeding Level

Estimates of genetic differentiation among populations from microsatellites ( ; ) indicated moderate but highly significant degree of differentiation among populations of Milicia excelsa in Benin. Since the species of concern is a tropical cross pollinated species with bat and bird dispersal-seeds [34] and some gravity-dispersed seed, this result may fit with the general observation that woody, perennial and outbreeding species maintain most of their variation within populations Hamrick et al., 1992 in [35]; [7]. Other studies that have assessed genetic variation in natural populations of tropical tree species using microsatellites exhibited much lower levels of genetic differentiation than found here (0.031 for Symphonia globulifera in Southern Costa Rica, [1]; 0.041 for Carapa guianebsis in Costa Rica, [4]; 0.032 for Swietenia humilis in Honduras, [8]; 0.097 for Swietenia macrophylla in the Brazilian Amazon, Lemes et al. [5]; and 0.026 for Vitellaria paradoxa in Mali, [33]). However, some authors have suggested that population differentiation is more accurately estimated by , because this measure better accounts for the high mutation rate of microsatellite markers while tends to underestimate population differentiation [9, 14, 36] Lemes et al. [5]. On the other hand, different simulation studies [37] have shown that most of the microsatellite markers formed by short or imperfect repeat motifs tend to fit more closely the IAM than the SMM.

Inbreeding coefficient ( ) indicated that alleles within populations were not united at random and that mating between close relatives may play an important role in determining the genetic structure of these populations. The estimated gene flow ( ) was also low reflecting a reduced gene flow among populations. The first possible explanation of such a situation is the occurrence of null alleles, which failed to amplify because of mutations in the flanking primer sequence (Lemes et al. [5]). But this seems to be an unlikely explanation because amplification failures were rare in this study. The second explanation which has been experienced in Swietenia macrophylla, a tropical outcrossing species is that, assortive mating caused by spatial clustering or coincidence in flowering time among related groups of trees, can lead to inbreeding and homozygote excess (Lemes et al. [5]). Such explanation seems to be more favorable in the case of Milicia excelsa in Benin. Indeed, Milicia excelsa is threatened throughout its range in Benin as a result of overexploitation and habitat destruction, which have clearly reduced local population size and led many populations to local extinction. The species’ populations currently occupy fragmented habitats such as sacred groves which are small patches of the original vegetation protected by traditional ethnobotany or religious practices, and farmlands where trees are protected by farmers through traditional agroforestry systems [38]. The spatial clustering of M. excelsa trees may have induced an increase local inbreeding level in progeny cohorts. Then the great majority of dispersal occurs between nearest-neighbour populations with only a few individuals moving longer distance. Other authors have also observed positive values and attributed it to populations’ substructure and inbreeding [33, 39].

5.2. Spatial Genetic Structure and Conservation Implications

Our data suggested that M. excelsa populations in Benin are differentiated by a process of “isolation by distance”. This hypothesis is supported by a significant Mantel test between a matrix of values and a matrix of geographical distances. Spatial analysis by the model-based clustering method also indicated a geographical structure and clustered populations into four groups which correspond to southern west, southern east (including one population from the middle northern east), northern west, and central populations. Several studies dealing with spatial genetic structure reported either a lack or a weak geographical heterogeneity in various tropical or temperate tree species which, they argued by the limited seed dispersal but extensive gene flow by pollen [4042] or by extensive gene flow by pollen, wide seed dispersal, self incompatibility and dispersal agent [43, 44]. These authors made the general statement that for woody insect-pollinated species with seed widely and independently dispersed by birds, at most weak genetic structure will result. In the current study, the distribution of gene pool 2 and gene pool 4 in the south and the predominance of gene pool 1 in the center, suggested an impact of limited seed dispersal. Different factors could explain this genetic structure. For instance, fragmentation and overexploitation may be a valuable explanation of this situation. As the remaining M. excelsa trees in these regions are grouped into small patches in the built-up village areas (fetish worship sites, open places, near-by-villages farms, [38]), dispersal agent activities and hence seed dispersal may be reduced by human disturbance. Bats that are the main pollinators of this species may tend to forage inside patches or between nearest-neighbor clumps. Hence fragmentation may cause genetic isolation by keeping these mammals inside fragments [9, 45], restricting gene flow among populations.

In addition, fruiting and seed dissemination of M. excelsa occurs during the dry season and germination rate is low and decreases rapidly. Therefore bulk part of the seed may loose its viability before it is transported far away by water flow during the rainy season. The mixed composition observed in Niaouli population harbouring mainly gene pools 2, 3 and 4 could be understood by the fact that this population is composed in majority of planted trees of unknown seed origin. We suspect that seeds used to establish plantations may have come from genetically distinct localities.

The information gained on the levels and distribution of microsatellite variation in Milicia excelsa could be used to suggest appropriate management strategies. M. excelsa has long suffered from habitat degradation caused by selective logging and land use pressure which is now inducing inbreeding in most of the studied populations and may compromise the species genetic diversity by genetic drift. Traditional strategies of conservation were helpful to protect some mother trees and to assure populations renewal but these practices are not efficient for genetic resources conservation. The high level of genetic variation found within populations combined with the low density of adults trees indicate that in situ conservation strategies should be designed to preserve large areas to minimize the loss of diversity due to genetic drift. The isolation by distance observed among populations suggests that in situ reserves should be distributed evenly across the species range in order to maximally conserve the regional diversity and integrity. Therefore we suggest that a key strategy for Milicia excelsa conservation should be the enhanced management of habitats to allow connections between populations. Agroforestry systems, if they are well managed, could contribute to maintain this connectivity across the landscape, hence acting as biological corridors between populations [46]. The species should be rehabilitated by transplanting seedlings into some preferential habitats such as Pobe, Sakete, Aplahoue and Bassila semi deciduous forests which should have high conservation priority as well as other protected reserves in savannah regions. These reserves could also serve as corridor among the existing patches. However, as mentioned by Storfer [47] and Gao [48], care should be taken to avoid mixing seeds or live plants collected from distinct populations to prevent changes in the genetic composition of the populations, producing a decreased fitness through outbreeding depression and disruption of locally adapted gene combinations.

Authors’ Information

This study is part of Christine Ouinsavi’s Ph.D. thesis on the silviculture, ecological structure, and conservation genetics of Milicia exclesa. Nestor Sokpon is a Professor interested in forest ecology and conservation of tropical tree species. Damase Khasa is a Professor interested in the molecular ecology of trees and their microsymbionts.

Acknowledgments

The authors are grateful to A. Gagné and S. Senneville (CEF, Université Laval, QC) for their technical assistance during the laboratory work. They also thank Dr. Jean-Luc Jany (Department of Plant Science, University of Cornell) for his helpful discussions on microsatellite work. This research was financially supported by Agence Universitaire de la Francophonie (AUF) and Bureau des Bourses et de l'Aide Financière (BBAF, Université Laval, QC) Fellowships to C. Ouinsavi, a competitive funds of Système National de la Recherche Agricole (SNRA, Benin) to N. Sokpon, and the National Sciences and Engineering Research Council (NSERC) of Canada discovery grant to D. P. Khasa.