Genetic and morphological diversity in populations of Annona senegalensis Pers. occurring in Western (Benin) and Southern (Mozambique) Africa

Study area

The study was carried out in Niassa Special Reserve (NSR), Mozambique (Southern Africa), and in the Sudanian zone, Benin (Western Africa), two locations where the species is best known and used. NSR is located in Northern Mozambique approximately between latitudes 12°8′40′N and 12°22′40′N; and longitudes 37°21′00′E and 37°45′00′E (Fig. 1). It covers approximately 42,000 km² and has been described as the largest protected area of Mozambique and the third largest in Africa (Ryan et al., 2016; Mbanze et al., 2019). Seventy-two percent of the total area of NSR is covered by dry Zambezian Miombo woodlands that are dominated by Brachystegia spiciformis Benth, Brachystegia boehmii Taub, and Julbernardia globiflora Benth. (White, 1983). The climate is tropical sub-humid, with a dry and relatively hot period between May and October. The annual rainfall is on average 900 mm per year increasing from the East (800 mm) to the West (1,200 mm). Temperature ranges between 20 and 30 °C (Allan et al., 2017). About 60,000 people are living inside the reserve and are concentrated around the two main villages of Mecula (Moz_MEC) in the East and Mavago (Moz_MAV) in the West, and along the main road (NCP, 2017; SRN, 2008; Mbanze et al., 2021). Slash-and-burn agriculture is the main livelihood activity of the population (Cunliffe et al., 2009).

The Sudanian zone is located in Northern Benin between latitudes 9°45′N and 12°25′N and longitudes 0°45′E and 3°55′E (Fig. 1) and is characterized by a tropical dry climate with two seasons (rainy and dry). The mean annual rainfall in this zone is often below 1,000 mm and the temperature is on average 27.5 °C (Gnanglè et al., 2011). The vegetation is composed of dry forests, woodlands, savannahs, and riparian forests. Common tree species in the area include Isoberlinia spp., Combretum spp., Acacia spp., Hyparrhenia spp., Loudetia spp., and Andropogon spp. (Gnanglè et al., 2012). North Borgou (Ben_BGN) and Mekrou pendjari (Ben_MPE) are the two main phytogeographical districts of the Sudanian zone of Benin. People living in the Sudanian zone of Benin are mainly farmers.

Sampling and data collection

A total of 154 individuals of Annona senegalensis from the two geographical regions (Fig. 1; Table 1) were analyzed for genetic diversity and population structure. Due to the unavailability of trees bearing mature fruits in some populations, morphological data focused on a total of 147 individuals. In each region, two sub-populations were selected and within each, and leaves and fruits were collected along a linear transect of 30 km, with a minimum distance of 100 m to 10 km to avoid sampling siblings. Twenty-seven to sixty individuals were sampled within each sub-population and used for genetic and morphological analysis (Table 1). Samples from all individuals were brought to the laboratory for morphological analysis. For the genetic analysis, fresh leaves were kept in silica gel while in the field and stored at −80 °C once in the laboratory until DNA extraction.

DNA extraction and nSSR amplification

Genomic DNA was isolated with the InnuSPEED Plant DNA Kit (Innuscreen GmbH, Analytik, Jena, Germany) according to the procedures described by the manufacturer, using 50 mg of ground leaves. The yield and purity of each sample was determined by spectrophotometry at 230, 260, and 280 nm (Nanodrop 2000, Thermo Fisher Scientific, Waltham, MA, USA), and complemented by agarose gel electrophoresis (Gomes et al., 2020). Based on the consistency of the polymorphic amplifications, ten microsatellite markers were used for genetic diversity analysis: LMCH4, LMCH6, LMCH11 (Escribano, Viruel & Hormaza, 2004), LMCH29, LMCH43, LMCH48, LMCH78, LMCH79, LMCH119, and LMCH122 (Escribano, Viruel & Hormaza, 2008).

PCR reactions were performed under the following conditions: 94 °C for 1 min; 94 °C for 30 s, 55 °C for 30 s (45 °C for LMCH29), 72 °C for 1 min (35 cycles); and 72 °C for 5 min. Each reaction was performed in a final volume of 15 µL containing 100 ng of genomic DNA, 0.4 µM each primer (Table 2), 1.25U MyTaq DNA polymerase and 1X MyTaq Reaction Buffer (Meridian Bioscience, Cincinnati, OH, USA). Forward primers were labeled with a fluorescent dye at the 5′-end. PCR products were separated by capillary electrophoresis on a CEQ™ 8000 capillary DNA analysis system (Beckman Coulter, Fullerton, CA, USA) and allele sizes were determined with GeneMapper 3.2 (Applied Biosystems, Waltham, MA, USA).

Data on morphological traits

Data were collected in Mozambique from January to April 2021 and in Benin from June to September 2021, after fruiting. Six morphological descriptors were measured on plants, namely total height, bole height, crown height, trunk diameter at breast height, crown diameter, and crown shape. The bole height is the height from the ground to the first big branch and the crown height is the difference between total height (m) and bole height (m). The crown shape was derived from the ratio of crown height over crown diameter. To determine the crown diameter (m), four radii were measured from the projection of the crown on the ground (Glèlè Kakaï et al., 2011; Hounkpèvi et al., 2016). At least 40 leaves and 40 ripened fruits were collected per individual. Seven morphological fruit descriptors were measured: fruit length (mm), fruit width (mm), fruit dry weight (g), number of seeds per fruit, seeds weight (g), pulp dry mass (g), fruit shape and the ratio fruit length to fruit width (Hounkpèvi et al., 2016; Lawin et al., 2021). Six quantitative descriptors were measured on leaves, including leaf length (cm), leaf width (cm), limb length (cm), petiole length (cm), leaf dry weight (g), and the ratio of leaf length to petiole length (Sun et al., 2020; Mollick et al., 2021). Fruits and leaves were further oven-dried at 105 °C until constant weight for the determinations of fruit dry weight, pulp mass, seeds dry weight, and leaf dry weight. After measuring the fruit dry weight, each fruit was split manually and the seeds were separated from the pulp. The number of seeds per fruit was then counted, and the seeds weight and pulp mass were weighed. Weights were measured using a 0.01 g precision scale while a centimeter ruler and a digital caliper with a 0.01 mm level of precision were used for all other measurements (Table 3).

Bioclimatic data

Using the GPS coordinates of each individual in QGIS 3.16.2 (QGIS Development Team, 2021), bioclimatic data was extracted from the CHELSA (Climatologies at High resolution for the Earth’s Land Surface Areas) database, considering the last data available over 30 years (1979-2013).

Genetic diversity, population structure, and differentiation

For each geographical area and sub-population, genetic diversity was assessed by calculating the total number of alleles (Ta), mean number of alleles per locus (Na), Shannon’s information index (H), mean expected heterozygosity (He), mean observed heterozygosity (Ho), inbreeding coefficient (FIS), and Polymorphism Information Content (PIC) using GenAlEx 6.51 (Peakall & Smouse, 2012). The Bayesian program STRUCTURE v.2.3.4 (Pritchard et al., 2000) was used to test whether any discrete genetic structure existed among samples. The analysis was performed assuming 1 to 10 genetic clusters (K) with ten replications per K. Models were run assuming ancestral admixture and correlated allele frequencies using run lengths of 300,000 interactions for each K after 50,000 burn-in steps. The optimum K was determined using STRUCTURE HARVESTER (Earl & VonHoldt, 2012), which identifies the optimal K based on both the posterior probability of the data for a given K and the ΔK (Evanno, Regnaut & Goudet, 2005). The results of the replicates at the best-fit K identified by STRUCTURE were then post-processed using CLUMPP 1.1.2 (Jakobsson & Rosenberg, 2007). A Principal Coordinates Analysis (PCoA) was also constructed in GenAlEx 6.51 (Peakall & Smouse, 2012) to detect the genetic relatedness among individuals based on Nei’s genetic distance. Analysis of molecular variance (AMOVA) was performed to quantify the partitioning of genetic variance between the geographical regions, as well as between and within sub-populations that showed genetic differentiation in STRUCTURE and PCoA. Each AMOVA was run with 10,000 permutations at 0.95 significance levels in Arlequin 3.11 (Exoffier, Laval & Schneider, 2005). The relationships between population pairwise Nei’s genetic distances and linear geographical distances (isolation by distance) were examined with a Mantel test (Mantel, 1967) implemented in Arlequin 3.11 (Exoffier, Laval & Schneider, 2005) using the same permutation and significance levels.

Morphological diversity and structuring

The morphological traits of 40 leaves and 40 fruits were measured per individual, and the respective averages were used for statistical analysis. Individual data was recorded for the remaining traits, i.e., total height, bole height, crown height, trunk diameter at breast height, crown diameter, and crown shape. The mean, standard error, and coefficient of variation of each morphological trait were calculated by population and sub-population. The coefficient of variation (cv %) was used to assess the variability of each morphological trait, considering a cv < 25% an indicator of weak variability (Reza et al., 2017). A student t-test was first used to evaluate differences between Northern and Southern populations. Similarly, an analysis of variance was used to compare traits among the four sub-populations. The assumptions of normality and homoscedasticity required to run these tests were checked previously, using the Shapiro–Wilks test and the Levene test, respectively. When the violation of the assumption of normality was severe (p < 0.01), the corresponding non-parametric test (Mann–Whitney or Kruskal-Wallis) was applied. When the Analysis of Variance (ANOVA) indicated a significant difference, a SNK-test was applied as a multiple comparison test in the package AGRICOLAE (De Mendiburu & Agricolae, 2020) to separate means.

To assess the relationship between the morphological descriptors and the bioclimatic variables, a redundancy analysis (RDA) within the VEGAN package was carried out on the least square mean values of the morphological descriptors and bioclimatic variables (Table S, supplementary data). The RDA was first carried out separately for the Western, and the Southern sub-populations. Another RDA analysis was implemented with the merged populations. These RDA analyses were intended to assess whether the relative importance of the relationships between bioclimatic variables and morphological variation was similar for the two regions. All analyses were implemented in R statistical software version 4.1.2 (R Core Team, 2021).

Genetic diversity, structure, and differentiation

A total of 156 alleles were found among the 154 Annona senegalensis samples. The number of alleles varied from 27 in the Southern (Moz_MAV sub-population) to 55 in the Western (Ben_BGN sub-population) region (Table 4). The total number of alleles was significantly higher in the sub-populations sampled in the Western region than in the Southern region (F = 3.23, p =0.023). This pattern was also observed in the average number of alleles (F = 2.05, p = 0.001), the Shannon Diversity Index (F = 1.04, p = 0.021), and the observed (F = 4.24, p = 0.019) and expected heterozygosity (F = 4.47, p = 0.024) (Table 4). The percentage of polymorphic loci was overall very high and showed the same pattern (F = 3.39, p = 0.025) i.e., higher in the Western than in the Southern population (Table 4). FIS showed negative values in all sampled populations (Table 4) suggesting a heterozygosity higher than expected under the Hardy-Weinberg assumption. FIS values were lower in the Western than in the Southern population (F = 1.29, p = 0.012; Table 4).

The Bayesian clustering program STRUCTURE found the highest LnP(D) and ΔK values for K = 2 differentiating the samples collected in Benin from the ones collected in Mozambique (Fig. 2). Nevertheless, STRUCTURE further revealed a secondary high LnP(D) and ΔK values at K = 4 differentiating the four sub-populations, Ben_BGN, Ben_MPE, Ben_MEC, and Ben_MAV into different genetic clusters (Fig. 2). Despite an overall high genetic integrity found in most samples, the results showed some signs of admixture between the genetic groups from Benin and Mozambique, although this admixture was negligible (Fig. 2). The same geographical pattern was retrieved by a principal coordinate analysis (PCoA) (Fig. 3). The first two coordinates of PCoA explained 35.9% of the total variation. Samples were spatially separated considering the two main geographic areas (Benin and Mozambique), but also by sub-populations following the K = 4 clustering result found in STRUCTURE (Fig. 2). The degree of spatial separation was lower for the two Mozambican sub-populations than for the ones from Benin (Fig. 3).

AMOVA revealed that a high proportion of genetic variation was attributable to significant differences between the two regions (69.1%) supported by high levels of genetic differentiation (FST = 0.305, p < 0.001). In addition, 21.3% of variation occurred among populations while the remaining was found within sub-populations. In addition, the Mantel test confirmed the existence of a significant positive correlation between Nei’s genetic distance and geographic distance for all pairwise sub-populations (r = 0.212, p < 0.001).

Morphological diversity and structure

The morphological traits of A. senegalensis varied significantly between Northern and Southern populations. Plants from Southern populations were significantly larger (DBH: 15.89 ± 2.10 cm) and taller (Total height: 5.72 ± 0.75 m) than those from Northern populations (DBH: 5.89 ± 0.62 cm; Total height: 2.56 ± 0.27 m) (Table 5). Irrespective of the population, the coefficient of variation (cv) was high (cv > 25%) for all traits. However, Northern populations had the highest cv values irrespective of the traits, except for the crown shape. A significant difference was reported among sub-populations. However, the sub-population Ben_BGN (North Benin) had a similar trunk diameter at breast height (DBH) to that of Moz_MEC (South Mozambique). Within Southern populations, sub-populations Moz_MEC and Moz_MAV had similar values for bole height. The cv value decreased from both populations (31.04% to 65.55%) to sub-populations (15.49% to 65.00%), but was still relatively high. The DBH, bole height, crown height, crown diameter, and crown shape were more diverse in Moz_MAV, while the total height highly varied in Ben_MPE (Table 5).

The morphological parameters of fruits and leaves varied significantly among populations (Table 6). Considering the four sub-populations, results showed that Moz_MAV had the highest value for fruit length, fruit shape, fruit dry weight, number of seeds, seeds weight, pulp dry mass, limb length, leaf length, and leaf width. Ben_BGN had the highest fruit width, while the highest value for petiole length and leaf weight was recorded in Ben_MPE. Fruits from Mozambique were found to be bigger than those from Benin. Some traits, like the ratio of leaf length to petiole length, showed similar values between populations (Table 6). Furthermore, Ben_BGN and Moz_MEC presented similar values for fruit length and number of seeds per fruit. Both regions and their respective sub-populations showed high cv regarding all traits, except fruit length, width, and fruit shape (Table 6). For fruit length, cv values varied from 15.53% to 15.96% in Western populations and from 14.53% to 18.76% in Southern populations. For fruit width, cv values varied from 11.68% to 13.43% in Western populations and from 12.46% to 14.58% in the Southern populations; and for fruit shape, cv values varied from 7.41% to 10.03% in Western populations and from 8.96% to 12.45% in Southern populations. The hierarchical clustering of the individuals based on their morphological traits resulted in four clusters (Fig. 4).

Overlap between genetic and morphological clusters

The chi-square test was performed to test the association between morphological clusters and genetic clusters. Results (Pearson chi-square = 209.771, DF = 9, p < 0.0001; Likelihood ratio chi-square = 195.358, DF = 9, p < 0.0001) suggested a significant association between genetic and morphological variation, and hence, an effect of genetic factors on the morphological variation i.e., the distribution of trees in genetic clusters is not independent of the morphological clusters. For instance, 86.67% of the individuals included in morphological cluster 1 corresponded to genetic cluster 1; and 70% of the individuals included in morphological cluster 3 correspond to the genetic cluster 3 (Table 7).

Influence of bioclimatic variables on the morphological variation

The redundancy analysis showed that there was a significant correlation between morphological traits and bioclimatic variables. Furthermore, this relationship varied in diverse ways according to the two regions. In all cases, only the first two axes were significant (p = 0.001, F = 12.489) and explained the extent to which variation in morphological traits is related to bioclimatic variables. In Western populations, the model considered nine out of the 19 bioclimatic variables (F = 7.7245, p = 0.001, adjusted R² = 0.404). The first axis (RDA1) explained 80.77% of the total variance and was a combination of mean diurnal air temperature range (chelsa_b_1), temperature seasonality (chelsa_b_3), mean daily maximum air temperature in the warmest month (chelsa_b_4), annual range of air temperature (chelsa_b_6), mean daily air temperature of the wettest quarter (chelsa_b_7), mean daily air temperature in the warmest quarter (chelsa_b_9), and mean monthly precipitation in the coldest quarter (chelsa_b_18). The second axis (RDA2) explained 8.56% of the total variation and combined mean annual air temperature (chelsa_bio), and annual precipitation (chelsa_b_11).

In Southern populations, the model considered only five out of the 19 bioclimatic variables (F = 5.3517, p = 0.001, adjusted R² = 0.279. The first axis (RDA1) explained 67.17% of the total variance and was a combination of mean annual air temperature (chelsa_bio), mean diurnal air temperature range (chelsa_b_1), temperature seasonality (chelsa_b_3), and mean daily maximum air temperature in the warmest month (chelsa_b_4). The second axis (RDA2) explained 21.74% of the total variation and only considered the annual range of air temperature (chelsa_b_6).

When merging the two populations, the model considered 11 out of the 19 bioclimatic variables (F = 12.489, p = 0.001, adjusted R² = 0.463). The first axis (RDA1) explained 83.44% of the total variance and was the combination of isothermality (chelsa_b_2); temperature seasonality (chelsa_b_3); mean daily minimum air temperature in the coldest month (chelsa_b_5); mean daily air temperature of the wettest quarter (chelsa_b_7); mean daily air temperature of the driest quarter (chelsa_b_8 ); annual precipitation (chelsa_b_11); amount of precipitation in the wettest month (chelsa_b_12); amount of precipitation in the driest month (chelsa_b_13); precipitation seasonality (chelsa_b_14); and mean monthly precipitation in the wettest quarter (chelsa_b_15). The second axis (RDA2) explained 6.19% of the variation and only considered the mean daily air temperature in the warmest quarter (chelsa_b_9). However, some of these bioclimatic variables such as mean annual air temperature, temperature seasonality, mean monthly precipitation in the coldest quarter, annual range of air temperature, and the mean monthly precipitation in the wettest quarter were not statistically significant in these different models (Table 8).

In the Western population, DBH, total height, bole height, crown height, crown diameter, crown shape, fruit length, fruit width, fruit shape, fruit dry weight, number of seeds per fruit, seeds weight, pulp dry mass, leaf width, petiole length, the ratio leaf length to petiole length, leaf dry weight and limb length were all loaded on RDA1, while only leaf length was loaded in RDA2 (Table 9). Based on the scores of morphological traits and bioclimatic variables on RDA axes (Tables 8, 9), leaf length was positively influenced by annual precipitation (chelsa_b_11). DBH, total height, bole height, crown height, crown diameter, fruit length, fruit width, fruit dry weight, number of seeds per fruit, seeds weight, and pulp dry mass were negatively influenced by mean diurnal air temperature range (chelsa_b_1), mean daily maximum air temperature in the warmest month (chelsa_b_4), annual range of air temperature (chelsa_b_6), and the mean daily air temperature of the wettest quarter (chelsa_b_7).

In Southern populations, all morphological parameters were loaded on RDA1 except bole height, crown diameter, crown shape, and the ratio leaf length to petiole length, which were loaded in RDA2. Crown diameter and the ratio leaf length to petiole length were negatively influenced by the mean annual air temperature (chelsa_bio), mean diurnal air temperature range (chelsa_b_1), temperature seasonality (chelsa_b_3), and the mean daily maximum air temperature in the warmest month (chelsa_b_4) (Tables 8, 9).

Considering the global set of sub-populations, morphological parameters were all loaded in RDA 1 except limb length, leaf width, the ratio leaf length to petiole length, and leaf weight. The DBH, total height, bole height, crown height, crown diameter, crown shape, fruit length, fruit width, fruit shape, fruit weight, number of seeds per fruit, seeds weight, pulp mass, and leaf length were positively influenced by the mean daily minimum air temperature in the coldest month (chelsa_b_5) and mean daily air temperature of the driest quarter (chelsa_b_8 ) (Tables 8, 9). The petiole length and the ratio leaf length to petiole length were negatively influenced by the isothermality (chelsa_b_2); temperature seasonality (chelsa_b_3); mean daily air temperature of the wettest quarter (chelsa_b_7); annual precipitation (chelsa_b_11) precipitation in the wettest month (chelsa_b_12); precipitation in the driest month (chelsa_b_13) precipitation seasonality (chelsa_b_14) mean monthly precipitation in the wettest quarter (chelsa_b_15) and the mean daily air temperature in the warmest quarter (chelsa_b_9), respectively (Tables 8, 9).