Genetic variation and cryptic lineage diversity of the Nigerian red-headed rock agama Agama agama associate with eco-geographic zones.

Nigeria is an Afrotropical region with considerable ecological heterogeneity and levels of biotic endemism. Among its vertebrate fauna, reptiles have broad distributions, thus, they constitute a compelling system for assessing the impact of ecological variation and geographic isolation on species diversification. The red-headed rock agama, Agama agama, lives in a wide range of habitats and, thus, it may show genetic structuring and diversification. Herein, we tested the hypothesis that ecology affects its genetic structure and population divergence. Bayesian inference phylogenetic analysis of a mitochondrial DNA (mtDNA) gene recovered four well-supported matrilines with strong evidence of genetic structuring consistent with eco-geographic regions. Genetic differences among populations based on the mtDNA also correlated with geographic distance. The ecological niche model for the matrilines had a good fit and robust performance. Population divergence along the environmental axes was associated with climatic conditions, and temperature ranked highest among all environmental variables for forest specialists, while precipitation ranked highest for the forest/derived savanna, and savanna specialists. Our results cannot reject the hypothesis that niche conservatism promotes geographic isolation of the western populations of Nigerian A. agama. Thus, ecological gradients and geographic isolation impact the genetic structure and population divergence of the lizards. This species might be facing threats due to recent habitat fragmentation, especially in western Nigeria. Conservation actions appear necessary.

Sub-Saharan Africa (SSA) has a high level of species diversity and endemism. However, the evolutionary mechanisms governing patterns of diversification in this region have continued to cause debate and controversy (Moritz et al. 2000; Hill and Hill 2001). Two main hypotheses have been advanced: (1) geographic isolation between populations restricted to rainforest patches (due to the Pleistocene forest refugia, montane, and riverine barrier model) driving allopatric speciation, and (2) ecological gradients associated with ecological speciation in parapatry (Endler 1977; Schluter 2000, 2001, 2009; Kirschel et al. 2011; Portik et al. 2017). While controversy remains over which model of diversification prevails for SSA taxa, both drivers may play varying roles in speciation and population divergence across a variety of Afrotropical taxa (Portik et al. 2017).

Pleistocene climatic fluctuations are important to several models of allopatric speciation in Afrotropical taxa (Leaché et al. 2016). The savanna and tropical forest habitats of SSA have a dynamic and prolonged history of expansion and contraction that intensified during the Pleistocene resulting in the broad expansion of savanna (Jacobs 2004; Dupont 2011; Leaché et al. 2016). These climatic oscillations generally drove the divergence of populations to produce new lineages, hence, affecting the spatio-temporal diversification patterns of vertebrate taxa (Portik et al. 2017). Accordingly, the oscillations markedly affected the geographic distributions of populations and the genomes of Afrotropical species.

The roles ecological gradients play in species diversification and population divergence in Afrotropical taxa are recognized increasingly. SSA contains several ecotones, and, thus, ecological factors potentially could influence evolutionary divergence and speciation (Smith et al. 1997, 2001, 2005; Freedman et al. 2010; Portik et al. 2017). The ecological gradient hypothesis, which posits that ecological speciation occurs as a result of divergent selection between populations in different environments, may explain patterns of some present-day diversity and reproductive isolation (Endler 1977; Schluter 2000, 2001, 2009; Kirschel et al. 2011; Portik et al. 2017). Accordingly, adaptation to different ecological niches can limit gene flow between populations and lineages in parapatry and in secondary contact after initial spatial separation (Rundle and Nosil 2005; Nosil and Feder 2009; Schluter 2009). Alternatively, divergent selection may also arise between sympatric populations occupying separate niches within a single geographic area (Rundle and Nosil 2005). Thus, ecological speciation can take place under any geographic arrangement (allopatric, parapatric, or sympatric) as long as divergent natural selection results in reproductive isolation. Accordingly, divergent selection can shape organismal phenotypic evolution and functional diversification and set the distributional limits for populations and species (Sexton et al. 2009; Arbour and López-Fernández 2014; Bossu and Near 2015). Evolutionary divergence due to ecological variation has been reported in vertebrates including non-avian reptiles (Ogden and Thorpe 2002; Thorpe et al. 2004, 2010, 2012; Freedman et al. 2010; Leaché et al. 2016) and birds (Smith et al. 2005; Kirschel et al. 2011).

Ecology is a central driver of patterns of evolutionary divergence and speciation (Couvreur et al. 2011). However, ecology plays contrasting roles in speciation models. For instance, conservation of ecological niches (niche conservatism) occurs when geographic isolation is driven by Pleistocene forest refugia, montane, and riverine barriers. In this case, taxa restricted to rainforest patches are unable to adapt to changing environmental conditions. Thus, niche conservatism could promote allopatric divergence in fragmented habitats by limiting adaptation to new environments when populations maintain an ancestral niche. In ecological speciation, selection to alternative ecological niches (divergent selection) would be a pre-requisite. In ecological speciation, variables comprising species’ niches would be expected to deviate less from a random geographic distribution over a phylogeny than in the case of geographic isolation arising due to the Pleistocene forest refugia, montane, and riverine barriers. Therefore, niche divergence may lead to lineage formation when populations adapt to new environments (Wiens 2004; Wiens and Graham 2005). In sum, niche conservatism and divergence in geographic isolation imply ecology is similar or the same between focal populations, while ecological speciation requires ecology to be different. However, both geographic isolation and local adaptation to ecological gradients in ecologically distinct habitats may be important drivers of evolutionary diversification in SSA.

An understanding of how ecological variation influences processes of divergence and eventually speciation will yield insights into the patterns of vertebrate diversification in SSA. Unfortunately, few such studies exist from this region. Examining the role of ecology in explaining the present-day species diversity requires the integration of genogeographic and ecological niche modeling (ENM) approaches on widely distributed taxa occurring in regions with high levels of ecological differentiation and species endemism.

Nigeria is an Afrotropical region that exhibits considerable ecological heterogeneity (Iloeje 2001) and, thus, has an ideal biota for testing hypotheses related to divergence owing to environmental adaptation. Several bioregions based on topography, climate, precipitation, vegetation, and soils occur there (Figure 1;Iloeje 2001). Nigerian ecological zones group broadly into forest and savanna. Humid tropical forests of the west and south have a longer period of rain accompanied by lower temperature, while savannas of the north have a shorter rainy season with relatively higher temperature (Aregheore 2009). In addition, clearing of the forests for agricultural purposes has created an expense of grassland savanna with few trees in some western areas. These ecological zones may influence patterns of diversity and distribution (Iloeje 2001).

Figure 1.

Map of Nigeria showing the different ecological bioregions; re-modified from Iloeje (2001).

Reptiles are well suited for assessing the impact of environment and geography on lineage formation and diversification because of their low vagility and strong responses to environmental factors (Fontanella et al. 2008). In Nigeria, some animals, for instance lizards, have wide ranging distributions spanning numerous ecological zones, yet others do not. Hence, testing of biogeographic hypotheses using such taxa can provide insights into how ecological variation shapes species’ distributions, genetic patterns, and evolutionary history.

Eight different species of Agama lizards have been described for Nigeria: A. agama, A. picticauda, A. lebretoni, A. parafricana, A. doriae, A. gracilimembris, A. paragama, and A. sankaranica (Nneji et al. 2018). Of these, the red-headed rock agama, A. agama, is the most widely distributed species (Nneji et al. 2018). Apart from Nigeria, this species occurs in most SSA countries including Benin, Burkina Faso, Cameroon, Cape Verde, Chad, Gabon, Ghana, Guinea, Guinea-Bissau, Liberia, Mali, Mauritania, Senegal, and Togo (Mediannikov et al. 2012). Sexual size dimorphism, dichromatism, and adult male breeding coloration are conspicuous traits of this species. Diagnosis often includes having a white underside, brown hind limbs, and a tail with a light stripe down the middle. The tail-stripe typically consists of about 6–7 dark patches along its side. Agama agama lives in social groups including a lead male, and few numbers of females and subordinate males. Lead males are usually aggressive and fight with other males to maintain their territories. Seasonal migration to breeding grounds or hibernation sites is common. However, populations of this species do not make extended movements of more than 200 km in any time of the year.

Although A. agama is well-suited to arid conditions, it inhabits a variety of habitats, including grassland, tropical forests, and urban settings (Gonçalves et al. 2012; Mediannikov et al. 2012; Nneji et al. 2018). The widespread distribution of this species suggests the presence of genetic diversity (Nneji et al. 2018), making it an intriguing species for study. Previous studies (e.g., Gonçalves et al., 2012; Leaché et al. 2016) associated the population histories and divergences in A. agama across SSA with habitat shifts due to climatic fluctuations during the Pleistocene. Studies of Leaché et al. (2014) recovered shallow levels of historical divergence (≤5 Ma) and weak relationships among populations within the West Africa A. agama complex, suggesting cryptic species diversity. Therefore, further studies are needed to investigate whether or not there is an ecological basis for local adaptation and speciation.

Herein, we test the hypothesis that ecology affected the current geographic patterns of distribution and genetic diversity of Nigerian A. agama using 2 approaches. First, we use genetic analyses [mitochondrial and nuclear DNA (nuDNA) sequences] to examine genetic structuring and diversification. Second, we use ENM based on climatic and geographic variables to investigate how environmental conditions may have shaped present-day genetic diversity. We also test the roles niche conservatism and divergence (i.e., ecology) played in lineage diversification. Specifically, we address 3 questions for Nigerian A. agama: (i) Does genetic structuring coincide with the geographic and ecological regions? (ii) Which abiotic factors may drive genetic diversification? (iii) Do lineages have similar or different ecological niches? Our results explore the role ecology plays in driving patterns of genetic variation and speciation.

Materials and Methods

Taxon sampling

Sampling procedures followed protocols approved by the Kunming Institute of Zoology Animal Care and Ethics Committee. Our sampling of 230 individuals of A. agama from 93 localities covered all ecological zones in Nigeria. For the genetic analyses, one individual each of closely related A. lebretoni and A. parafricana were collected for outgroup analyses. On capture, a small tail biopsy (about 1.0 cm) was taken from each lizard and preserved in 95% ethanol. Subsequently, some animals were released at the site of collection. Fifty-eight voucher specimens were euthanized and tissues (liver/muscle) were sampled and stored in 95% ethanol. These lizards were fixed using 4% formalin and then transferred to 70% ethanol. Vouchers were deposited in the Department of Zoology, University of Ibadan and Gashaka Gumti National Park, Nigeria. Ethanol-preserved tissues were stored at −80°C in the Animal Genetic Resource Bank, Kunming Institute of Zoology, Chinese Academy of Sciences, China. Locality information and GenBank accession numbers were summarized in Supplmentary Table S1.

Molecular techniques

We extracted total genomic DNA from frozen tissues using standard protocols (Sambrook et al. 1989). Partial sequences of mitochondrial DNA (mtDNA) cytochrome c oxidase subunit I (COI; Che et al. 2011) were sequenced for all individuals. We also amplified and sequenced nuDNA markers, including oocyte maturation factor MOS (CMOS) (Saint et al. 1998) and RNA fingerprint protein 35 (R35) (Leaché 2009). The nuDNA sequences were acquired from a subset of individuals (n = 80) selected from the matrilineal genealogy (see the “Results” section). We also sequenced mitochondrial and nuDNA markers of the closely related A. lebretoni and A. parafricana. Primer sets used for amplification via the polymerase chain reaction (PCR) and Sanger sequencing were shown in Supplementary Table S2. Amplification was performed in 25 µL volume reaction that contained 1.5–2.5 µL of genomic DNA, 18.5 µL HOH, 2.5 µL of PCR buffer, 1 µL of dNTPs, 1 µL of each of the primer (10 pm/µL), and 0.30 µL of r/l-Taq polymerase. PCR cycle profiles were as follows: an initial denaturation for 5 min at 94°C, followed by 35 cycles of 1 min at 94°C, annealing for 45 s at 46–61°C, extension for 1 min at 72°C, and a final extension for 10 min at 72°C. The resulting PCR products were visualized on 1.2% agarose gel containing gel red. Purified PCR products were directly sequenced with an automated DNA analyzer (ABI 3730 xl) in both forward and reverse directions.

Sequence alignment and partitioning

The nucleotide sequences were viewed and confirmed by eye using SeqManTMII (DNASTAR Lasergene 7). Alignment used ClustalW (Tamura et al. 2013) with default parameters in MEGA 6.0 (Tamura et al. 2013). The aligned sequences were translated into amino acids to check for premature stop codons and to confirm that the open-reading frame was maintained in the protein-coding loci. All sequences were submitted to BLAST searching in NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi) to confirm identity. Nuclear gene sequences containing more than one ambiguous site were resolved using PHASE (Stephens et al. 2001) as implemented in DNASP v5.10 (Librado and Rozas 2009) with default parameters. Haplotypes for mtDNA and alleles of the phased nuDNA were generated using DNASP. Best-fit partitioning schemes and substitution models were identified for each gene separately using PartitionFinder v1.1.0 (Lanfear et al. 2012). All PartitionFinder analyses used the following searching criteria: branch lengths =linked; search =greedy, models =mrbayes, model_selection =BIC; datablocks partitioned for 1st, 2nd, and 3rd codon positions for each of the genes; and search =all. The PartitionFinder selected partitioning of mtDNA data into Codon 1+2, and 3, and nuclear data unpartitioned by codon, but partitioned by region/locus (i.e., different partition for CMOS and R35).

MtDNA genogeography and species tree estimation

To answer our first question as to whether or not the genetic structuring of Nigerian A. agama coincided with the geographic and ecological regions, we estimated the mtDNA genogeography and species tree using our molecular data. Our matrilineal genealogy was inferred from 31 unique haplotypes using Bayesian inference (BI) as implemented in MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003). The analysis used partitions and models identified with PartitionFinder, and was run for 20 million generations with sampling every 1,000th generation. Two independent runs with four MCMC chains were performed. We excluded the first 25% of trees as burn-in before the log-likelihood scores stabilized. A 50% majority rule consensus of the sampled trees was constructed and visualized using FigTree 1.4.2 (Rambaut 2012). We considered bootstrap values of Bayesian Posterior Probabilities (PPs) ≥0.95 as being strongly supported (Hillis and Huelsenbeck 1992).

A species tree was estimated under a coalescent model using BEAST 1.7 (Drummond et al. 2012). The analysis contained 80 samples belonging to the 4 matrilines identified by our matrilineal genealogy (“see the Results” section) plus the outgroup species. Our mtDNA gene tree was used to assign individuals to groups. The site models, clock models, and gene-trees were unlinked across the three loci. We employed a Yule tree prior and strict clock model that allowed every branch in the tree to evolve according to the same evolutionary rate. Analyses were run 100 million generations, with parameters sampled every 10,000 generations. We checked convergence and stability of the parameters using TRACER 1.6 (Rambaut et al. 2014). The species tree was summarized after discarding the first 25% of trees as burn-in.

Population genetic structure and demographic analyses

Because our first question sought to unravel genetic structuring, we further used additional statistics to investigate patterns of population structure. We used uncorrected pairwise distances to estimate the evolutionary sequence divergences between and within matrilines identified by our mtDNA genealogy (see the “Results” section) using MEGA 6.0 (Kumar et al. 2008; Tamura et al. 2013). We calculated the pairwise Fst values between and within matrilines using ARLEQUIN 3.5 (Excoffier and Lischer 2010). Further, we estimated diversity indices by computing number of haplotypes (H), haplotype diversity (Hd), nucleotide diversity (π), and mean number of pairwise differences (k) for each locus and matriline in DNASP (Librado and Rozas 2009). We tested the null hypothesis of Hardy–Weinberg expectations (HWEs) using the nuDNA data set. Chi-square tests for significant deviations from HWE were conducted with ARLEQUIN.

To test if genetic differentiation between populations was owed to isolation-by-distance (IBD), which required a significant correlation between geographic and genetic distance (Fst), we performed a Mantel test because of the absence of established alternative methods that can detect the effect of non-independent geographical variables on gene flow (Guillot and Rousset 2013) and also, given that partial Mantel test has been recently criticized because of high type-I error rate and low resolution power (Raufaste and Rousset 2001; Harmon and Glor 2010; Legendre and Fortin 2010; Guillot and Rousset 2013). A Mantel test for matrix correlation between pairwise genetic distance and geographic distance was performed using the IBD Web Service (IBDWS) (Jensen et al. 2005). We assessed the significance using 10,000 matrix permutations.

Further, because of the obvious mitochondrial population structure identified by our mtDNA genealogy (see the “Results” section), we investigated evidence of population genetic structuring using Analysis of Molecular Variance (AMOVA; Excoffier et al. 1992) as integrated into ARLEQUIN. The grouping of populations was based on our matrilineal genealogy. Permutation tests for hierarchical AMOVA were performed at 3 levels: among groups, among populations within groups, and within populations, and populations grouped based on the type of vegetation zones. To test population structuring based on ecological type, we formed 3 groups of populations: forest specialist, derived savanna/forest specialist, and savanna specialist. Further, we tested geographic population structure based on the geographic regions. To do this, we formed 3 groups: North/East, South, and West. Significance was assessed by 10,000 random permutations.

To identify divergent genetic groups that might correspond to potential cryptic species within populations of Nigerian A. agama, we performed a Bayesian population mixture analysis implemented in BAPS 5.3 (Corander and Tang 2007; Corander et al. 2008), following the recommendations of Corander et al. (2009). BAPS assigned individuals to distinct gene pools probabilistically, based upon multilocus genetic data. BAPS does not make a priori assumptions about the number of gene pools (k) but a fixed number could have been set. We conducted the genetic mixture analysis using 78 individuals that had been sequenced for COI, CMOS, and R35. We performed “clustering of individuals” analyses using independent loci model in 2 ways. First, we forced BAPS to partition the individuals into 4 groups (k = 4), each reflecting a matriline. Then, we ran BAPS with maximal number of groups (k) set from 2 to 10 to determine the most probable number of distinct gene pools. For both searches we used 10 replicates.

Given the distinct matrilines identified by our mtDNA genealogy, we examined the demographic histories of each matriline using an array of statistics based on the COI data. We excluded matriline A (southern populations) as it was represented by a single haplotype. We conducted Tajima’s D (Tajima 1989), Fu and Li’s D and F (Fu and Li 1993) tests using ARLEQUIN. We assumed that if the population sizes had been stable across time, Tajima’s D and Fu and Li’s D would be near zero. Significantly positive values were assumed to result from recently experienced bottlenecks, and significantly negative values indicated recent population expansions (Tajima 1989; Fu and Li 1993). Because causation is difficult to ascertain when Tajima’s D deviates significantly from zero, we also computed Fs (Fu 1997). This statistic appeared to be particularly useful for detecting population expansions. A negative value of Fs was expected if a recent population expansion or genetic hitchhiking occurred and a positive value from a recent population bottleneck. Further, we examined for signals of population expansion by employing mismatch distributions (Rogers and Harpending 1992) as implemented in ARLEQUIN. A multimodal distribution was assumed to indicate demographic stability, while expansion generated a unimodal pattern (Slatkin and Hudson 1991). We compared observed distributions of nucleotide differences between pairs of haplotypes with those expected under spatial (Excoffier 2004) and demographic (Excoffier and Lischer 2010) expansion models by using the generalized least square approach. The sum of squared deviations (SSDs) were used as goodness-of-fit statistics for the observed and expected mismatch distributions, and the significance of the fit of the expansion model was tested, while the confidence intervals for the associated parameters estimates using 1,000 bootstrap replicates were examined using ARLEQUIN.

Ecological niche modeling

To answer our second question that sought to investigate how environmental conditions may have shaped patterns of present-day genetic diversity, and which abiotic factors may drive population genetic diversification, we employed ENM with Maxent 3.3 (Phillips et al. 2010) using 93 occurrence localities, 25% random points, and 19 bioclimatic variables were downloaded from the WorldClim 1.4 database at a spatial resolution of 2.5 m (Hijmans et al. 2005). Bioclimatic grids were clipped to an extent covering Nigeria. Replicate runs (n = 25) of MaxEnt using the “subsample” method were performed to ensure reliable results. Model performance was assessed using the area under the curve (AUC) of the receiver operating characteristic plot, with 25% of the localities randomly selected to test the model. AUC scores could have ranged between 0.5 (randomness) and 1 (exact match), with those above 0.90 being considered excellent (Swets 1988; Elith 2002). The MaxEnt jackknife analysis was used to evaluate the relative importance of the 19 bioclimatic variables employed, based on their gain values when used in isolation.

To test the roles niche conservatism and divergence played in lineage diversification and to further examine the extent of ecological differentiation between lineages, we performed the background similarity test. Specifically, we performed background similarity test on the two western matrilines identified by our matrilineal genealogy. To evaluate if potential ecological niches of both western lineages (matrilines B and C) and north/eastern lineage were more different from one another than expected based on the differences in the environment they occurred, we used Schoener’s D statistic (Warren et al. 2008) for niche overlap and the similarity statistic I (Schoener 1968) to quantify niche similarity. Both statistics could have ranged from 0 (no overlap) to 1 (identical niche models). Bootstrapping with 100 pseudoreplicates was used. Niche conservatism was expected if the observed niche overlap was significantly greater than the simulated overlap values while niche divergence was expected if the observed niche overlap was significantly less than the simulated overlap values.

Point-based analysis

To answer our third question—do lineages have similar or different ecological niches?—we performed a point-based analysis. The values for each of the 19 bioclimatic variables for all sampling sites were extracted using DIVA-GIS (Hijmans et al. 2012). The principal component analysis (PCA) was used to reduce the number of explanatory variables by creating new compound variables as PC1, PC2, PC3, and PC4. One-way MANOVA (Wilks’s λ) was conducted with resulting PC1, PC2, PC3, and PC4 as the dependent variable and the four matrilines as class variable to test if ecological niches differed significantly. To determine which dependent variables were responsible for significant MANOVA results, ANOVA tests were performed. Additionally, we used post hoc tests (i.e., Tukey’s HSD) to examine if the observed significant differences were between all matrilines or certain pairs. Analyses were performed using XLSTAT.

Results

Sequence information

The mtDNA (COI) dataset comprised 619 bp of aligned sites. Among these sites, 118 were polymorphic while 107 were potentially parsimony-informative. The gene consisted of 31 unique haplotypes, with Hd =0.94, nucleotide diversity (π = 0.460), and mean number of pairwise differences (k = 25.18). The two nuDNA datasets comprised 335 bp of aligned sites for CMOS and 632 bp for R35. Observed potentially informative sites were 2 (CMOS) and 17 (R35). The number of haplotypes (H) were 4 (CMOS) and 28 (R35). Overall, Hd, nucleotide diversity (π), and the mean number of pairwise differences (k) were as follows: CMOS (Hd =0.49, π = 0.014, k = 0.54); and R35 (Hd =0.85, π = 0.053, k = 1.886). Expected Hardy–Weinberg heterozygosity (HE) was higher in R35 (HE =0.85244) than in CMOS (HE =0.49096) (Supplementary Table S3). In both loci, observed heterozygosity was significantly less than expected values (P < 0.001).

MtDNA genogeography and species-tree

The BI analyses of the mtDNA haplotypes (Figure 2II) recovered four unambiguous matrilines (PP =1.00) whose distributions corresponded to geography and ecological (vegetation) zones (Figure 2I). Several well-supported sub-lineages (PP ≥0.95) showed strong concordance with ecology. These were as follows:

Figure 2.

(I) Map of Nigeria showing geographic distribution of matrilines of Agama agama. (II) Bayesian 50% majority-rule consensus tree of Nigerian A. agama inferred from mitochondrial COI sequences. Values above branches are Bayesian posterior probabilities (PP≥0.95); values below PP<0.95 not shown. Letters A–D above the branches indicate matriline and letters with numbers indicate sublineages. (III) Species-tree resulting from the *BEAST coalescent analysis based on mtDNA (COI) and nuDNA (CMOS and R35) sequences of Nigerian A. agama.

Matriline A (South): restricted to humid tropical rainforest of southern Nigeria.

Matriline B (West-1): consisted of 2 sub-lineages distributed across the humid tropical rainforests of Western Nigeria.

Matriline C (West-2): occurred in western Nigeria, sub-matrilines C1 and C2 were found within humid tropical forests, and C3 within the region historically referred to as derived savanna.

Matriline D (North/East): the widest distribution occurring within savannas of northern and eastern Nigeria and with well-supported sub-matrilines D1–D4.

The species-tree (Figure 2III) recovered lineages with relationships typical of our matrilineal genealogy, and particularly for lineages A and B, although with weak support (PP = 86). Also, the species-tree topology depicted lineages C and D as sister-groups (Figure 2III), but again with little support (PP =57), and the matrilineal genealogy failed to recover sister-group relationships for these lineages (Figure 2II).

Population genetic structure and historical demography

MtDNA-based estimates of the evolutionary divergence (uncorrected pairwise distance) and pairwise Fst values were shown in Table 1. Moderate mean genetic divergences (5.10–7.00%) occurred between matrilines (Table 1). The pairwise Fst values between matrilines ranged from 0.672% to 0.899%. The mtDNA divergence between the western matrilines B and C exhibited moderate genetic differences (5.50%) and pairwise Fst value (0.672%). The number of haplotypes (H) within matrilines (excluding A) ranged from 5 (matriline B) to 16 (matriline D) (Supplementary Table S4). Hd was high for all matrilines, except for B. Nucleotide diversity (π) and the mean number of pairwise differences also varied among matrilines (Supplementary Table S4).

Table 1.

MtDNA-based estimates of evolutionary divergence [uncorrected pairwise distances; (%) and population pairwise Fst values between matrilines of A. agama]

Group	Matriline A	Matriline B	Matriline C	Matriline D
Matriline A (South)	*	5.20	6.94	7.00
Matriline B (West 1)	0.899	*	5.50	5.10
Matriline C (West 2)	0.759	0.672	*	6.55
Matriline D (North/East)	0.865	0.783	0.726	*

Note: Top of the matrix is the uncorrected p-distances; below the matrix is the population pairwise Fst values.

Strong IBD was apparent for the entire set of 93 populations. The Mantel test was significant [correlation coefficient (r) =0.4793, P < 0.001, (Supplementary Figure S1). Thus, the null hypothesis of no correlation between mtDNA genetic and geographic distances was rejected. Further, our AMOVAs revealed significant variation among groups and high ΦST values (Supplementary Table S5), thus indicating high levels of mtDNA geographical structuring.

The initial analysis of genetic mixture by BAPs (k = 4) revealed four genetically differentiated eco-geographic groups (marginal likelihood =−1943.4157; Supplementary Figure S2). This included population clusters typical of our matrilineal genealogy (Supplementary Figure S2). With BAPS allowing the optimal value of k > 4, the analysis returned k = 6 (marginal likelihood =−1831.4551; Supplementary Figure S2) as the clear optimal value for the dataset. Generally, this conformed to eco-geographic differentiation with two groups partitioned as forest specialist (“Southern” and “Western” groups), two groups shared by individuals from derived savanna/forest region of western Nigeria, and the other two groups shared by populations from savanna region of northern/eastern Nigeria.

In terms of signatures of population expansion of the matrilines, overall, non-significant positive Tajima’s D, Fu’s and Li’s D and F statistics were obtained for matrilines B–D, except that matriline D had a non-significant negative value of Fu and Li’s D (Supplementary Figure S3). Non-significant positive Fs values ranged from 1.371 (matriline D) to 13.173 (matriline C) (Supplementary Figure S3). Mismatch distributions for all matrilines were multimodal (Supplementary Figure S3). However, the mismatch distribution test statistics (SSD and r) were relatively small and not significant (Supplementary Figure S3).

The ENM had excellent predictive power for occurrences. The four matrilines exhibited high sensitivity and specificity (>0.95 for all matrilines; Figure 3) and differed from the random prediction (AUC =0.5). The variables with the greatest contributions to the model for each lineage were as follows: Matriline A, South (minimum temperature of coldest month; 47.30%); matriline B, West 1 [temperature annual range (P5–P6); 42.70%]; matriline C, West 2 (precipitation of driest quarter; 40.10); and matriline D, North/East (precipitation of the warmest quarter; 21.90%) (Table 2). Thus, temperature variables were the highest ranked environmental variables for forest specialist while precipitation variables were the highest ranked variable forest/derived savanna, and savanna specialists (Table 2).

Figure 3.

Species distribution model for the matrilines of Nigerian A. agama based on current climate observations. Probability of occurrence is represented by different colors from low (blue) to high (red).

Table 2.

The estimate results of relative contributions of the environmental variables to the MaxEnt model for the matrilines of A. agama

	Percentage contribution
Environmental layer	Matriline A	Matriline B	Matriline C	Matriline D
Min temperature of coldest month	47.3	0	0.5	2.2
Temperature annual range	25.3	42.7	1.4	14.2
Precipitation of driest quarter	9.2	0	40.1	12.3
Isothermality	6.9	0	2.3	6.5
Precipitation of driest month	5.3	6.1	12.9	17.4
Precipitation of coldest quarter	4.8	0	12.7	3.7
Precipitation of wettest month	0.8	0	0	0
Annual precipitation	0.6	1.6	0	0
Mean temperature of driest month	0	0	1.7	0
Mean temperature of wettest quarter	0	0	0	1.0
Max temperature of warmest month	0	1.8	0	0
Temperature seasonality	0	11.2	24.7	0.4
Mean diurnal range: mean of monthly	0	7.1	0	0
Precipitation of warmest quarter	0	3.6	0	21.9
Precipitation of wettest quarter	0	0.9	0	0
Precipitation seasonality	0	21.7	0	0
Mean temperature of coldest quarter	0	0	2.5	1.2
Mean temperature of warmest quarter	0	0	0	13.2
Annual mean temperature	0	0	1.3	6.1

The background similarity test revealed a statistically significant similarity of abiotic niches between matriline C and the background of matriline B; values of Schoener’s D and similarity statistic I were significantly higher than the null distributions obtained from 100 randomizations (P < 0.01; Supplementary Table S6). The opposite relationship, matriline B versus the background of matriline C, was not significant (Supplementary Table S6). However, pairwise comparison of ecological niches showed that the geographically widely separated north/eastern and western matrilines differed strongly in their realized climatic niches (Supplementary Table S6). These results suggested that the climatic niches of matrilines B and C did not differ significantly, suggesting that niche conservatism plays a role at small geographic scales.

Point-based analyses

The PCA reduced the 19 bioclimatic variables to four PCs (Table 3). The first 4 PCs cumulatively explained more than the 96% of the variation, with each PC accounting for 60%, 18%, 10%, and 8% of the variation. PC1 explained variation from temperature variables (minimum temperature of the coldest month and isothermality) (Table 3) and PC2 was correlated with variables describing mean temperature of wettest quarter, annual temperature, and mean temperature of warmest quarter) (Table 3). PC3 loaded mainly for precipitation of the wettest quarter as well as precipitation of wettest month, and PC4 was dominated mainly by precipitation of the wettest month and precipitation of the driest month. The MANOVA test revealed significantly different ecological niches between the matrilines of Nigerian A. agama (Wilks’ lambda =0.027; F12,228 =55.58; P < 0.001; Supplementary Table S7), and ANOVA tests revealed that the association of the four PCs were significant (P < 0.001). Further, the post hoc Tukey’s test showed that significant (P < 0.05; Supplementary Table S7) pair-wise comparisons for certain pairs for PC1–3. Thus, our analyses showed that all matrilines were ecologically differentiated for at least one of the three PCs (PC1–3; Supplementary Table S7).

Table 3.

The PCA of the 19 environmental variables

	Name of environmental data	PC1	PC2	PC3	PC4
bio1	Annual mean temperature	0.430	0.897	−0.040	−0.042
bio2	Mean diurnal range: mean of monthly (max temp − min temp)	−0.957	−0.016	0.181	−0.097
bio3	Isothermality: (bio2/bio7)×100	0.929	−0.196	−0.150	0.037
bio4	Temperature seasonality	−0.887	0.101	−0.096	0.409
bio5	Max temperature of warmest month	−0.848	0.466	0.130	−0.109
bio6	Min temperature of coldest month	0.960	0.180	−0.178	−0.056
bio7	Temperature annual range (P5–P6)	−0.977	0.061	0.170	−0.005
bio8	Mean temperature of wettest quarter	0.222	0.902	−0.118	0.218
bio9	Mean temperature of driest quarter	0.820	0.451	−0.166	−0.251
bio10	Mean temperature of warmest quarter	−0.348	0.896	0.003	0.203
bio11	Mean temperature of coldest quarter	0.721	0.563	0.153	−0.348
bio12	Annual precipitation	0.898	−0.068	0.381	−0.014
bio13	Precipitation of wettest month	0.484	0.020	0.719	0.456
bio14	Precipitation of driest month	0.826	−0.033	−0.195	0.451
bio15	Precipitation seasonality (coefficient of variation)	−0.878	0.138	0.256	0.359
bio16	Precipitation of wettest quarter	0.561	0.054	0.797	0.186
bio17	Precipitation of driest quarter	0.869	−0.095	−0.244	0.393
bio18	Precipitation of warmest quarter	0.806	−0.215	−0.166	0.366
bio19	Precipitation of coldest quarter	0.709	−0.123	0.482	−0.448

Note: Max denotes the maximum value and min denotes the minimum value.

Discussion

This is the first study to explicitly use genetic data and ecological niche models to examine patterns of diversification of Nigerian A. agama. Analyses provide insights into the role ecological and geographic variation play in shaping the species’ diversification. Few studies have attempted to do this for other West African species. Hence, our results provide important insights into the biogeography of Nigeria.

Cryptic diversity and patterns of diversification

Substantial genetic differentiation and population structuring occurs within Nigerian A. agama. The four strongly supported matrilines have deep divergences (Figure 2II) and several factors may explain this.

Ecological variation appears to play a significant role in the genetic structuring of Nigerian A. agama. The lower Guinean forest, which also occurs in Nigeria, contains several ecotones that promote speciation through disruptive selection (Smith et al. 1997; Portik et al. 2017). Our results show that the matrilines are restricted ecologically to forests, forest/derived savanna, and savanna. ENMs also predict potential distributions with high confidence, indicating that abiotic environmental factors substantially influence the geographic ranges of the matrilines. Given that the main contributing parameters are “temperature variables” for the forest specialist and “precipitation variables” for the forest/derived savanna, and savanna specialists, this seems reasonable. Adaptation to local environmental conditions is a leading driver of evolution and speciation (Zhang et al. 2014). Adaptation of the matrilines to different local environmental conditions appears to have led to population genetic divergence. This implies differential abilities of populations to survive and reproduce in varying local environmental conditions. Thus, our analyses cannot reject the hypothesis that habitat differentiation together with the antecedent environmental variables can lead to speciation and diversification of animal in the Afrotropics (Smith and Bernatchez 2008). Our finding is congruent with other studies of Afrotropical taxa (Smith et al. 1997, 2001, 2005; Swart et al. 2009; Freedman et al. 2010). Similar results have also been reported in the rainforest skink Trachylepis affinis of SSA (Pinho et al. 2007; Freedman et al. 2010; Miraldo et al. 2011; Milá et al. 2013).

Geographic isolation over time seems a likely factor in the diversification of Nigerian A. agama. The IBD analysis obtained a significant correlation between mtDNA divergences and geographic distances. In addition, our hierarchical AMOVAs (Supplementary Table S5) revealed significant variation and high ΦST values (ΦST =0.771) when populations are grouped based on geographic location, thus indicating high levels of mtDNA geographical structuring. Quaternary climatic oscillations have shown to have profoundly affected the geographic occurrences and genetic structure of many African species (Portik et al. 2017). Plio–Pleistocene glaciations have been reported to influence the evolutionary history of many African animals such as amphibians (Portik et al. 2017), birds (Voelker et al. 2010), rodents (Nicolas et al. 2011), and primates (Telfer et al. 2003; Anthony et al. 2007). Previously, Leaché et al. (2014) dated divergence (≤5 Ma) within the West African A. agama complex. Hence, it seems likely that Plio–Pleistocene climatic fluctuations influenced the diversification and geographic occurrences of Nigerian A. agama. These events highlight the key role climatic fluctuations, and the associated expansions and retractions of forests, played in the evolution Africa’s faunal diversity (Portik et al. 2017).

Our analysis cannot test the hypothesis that human activities are a contributor to the diversification of Nigerian A. agama. However, matriline C possesses a signal of a recent population bottleneck or sharp reduction in population size. Likewise, Leaché et al. (2016) inferred reduction in range size in the West African A. agama complex during the Pleistocene. Individuals belonging to matriline C occur in the historically derived savanna region, which was originally the drier part of the high forest (Iloeje 2001). Changes in vegetation cover and associated climatic change appear to have influenced its diversity. Our analyses resolve a good fit to a multimodal distribution, which indicates a structured or stable population over time (Castoe et al. 2007). However, mismatch distribution test statistics (SSD and r) are quite small and not significant (Supplementary Figure S3). This pattern suggests an expanding population that has undergone some historical contraction (Castoe et al. 2007).

Analyses show that heterogeneous habitats and habitat-mediated climate differentiation probably acted as important physical barriers between regional populations of Nigerian A. agama, thus generating well-differentiated matrilines. The geographic distribution of the matrilines adds to the growing body of evidence that ecological variation and geographic isolation drive diversification in tropical species.

Niche conservatism versus divergence

Rather than finding that a single process of niche evolution played a dominant role in the diversification of Nigerian A. agama, our results show that different modes of niche–barrier interaction associate with lineage formation. Lizards have provided classical examples of evolutionary divergence frequently overcoming niche conservatism (Knouft et al. 2006; Losos 2008). Our results provide evidence that niche conservatism, rather than niche divergence, is the main driver of diversification in the two western matrilines. This is not unusual. Niche conservatism has been shown to influence the distribution of many taxa (Weins et al. 2006, 2009; Pyron and Burbrink 2009). Niche differentiation appears to be based on the dominant ecological feature of the local environment. Thus, while temperature exerts the strongest influence on matriline B, which inhabits the forest ecosystem of western Nigeria, precipitation appears most likely to affect matriline C, which occurs in the forest/derived savanna ecosystems of western Nigeria.

The mechanistic explanation for diversification within the western matrilines may relate with the inability of populations to adapt quickly to novel climatic conditions. Either stabilizing selection, gene flow, genetic constraints due to pleiotropy, and/or lack of appropriate genetic variation (Wiens and Graham 2005) may not have allowed western populations of Nigerian A. agama to adapt to Pleistocene conditions. However, we cannot rule out effects of niche divergence in promoting the diversification, especially for populations from southern and north/eastern Nigeria; it remains possible that matrilines from southern and north/eastern Nigeria exhibit divergence in environmental preferences on short evolutionary timescales. This possibility requires further exploration to determine how the two mechanisms drive diversification.

Implications for taxonomy and conservation

With regard to taxonomy, one possible interpretation of the genealogical pattern, particularly ecological and geographic divergence, is that A. agama is a complex of several cryptic species. Our result was consistent with the previous study of Leaché et al. (2009) on the West African A. agama group, which speculated that A. agama (sensu lato) contained multiple distinct lineages that co-occur throughout West Africa. Moreover, congruent with Leaché et al. (2016), the evolutionary divergence within A. agama may indicate speciation due to ecological variation and/or geographic isolation of populations. The geographic structure appears to result from long-term isolation and/or local adaptation. This awaits testing for gene flow by using large number of loci. Regardless, analyses clearly show that this species has multiple management and conservation units. Each lineage should be considered independently when setting up strategies for long-term conservation.

Authors’ Contributions

L.M.N., A.C.A., A.A.A.U., M.-S.P., R.W.M., and J.C. designed the project. L.M.N., A.C.A., A.O.O., O.C.O., Y.S., D.S., A.O.A., A.B.O., T.E.O., I.C.N., O.O., S.O.O., O.A.I., J.Y.U., T.A., A.D.U., H.A., O.O.I., J.Q.J., O.A.U., and A.A.A.U. aid in field survey and sample collection. L.M.N. performed the laboratory work and wrote the manuscript. L.M.N., A.C.A., and F.Y. analyzed and interpreted the data. A.C.A., A.A.A.U., M.-S.P., R.W.M., and J.C. edited the manuscript. All authors read and approved the last version of the manuscript.

Click here for additional data file.