Introduction

Mangrove forests are among the most dynamic habitats in tropical and subtropical intertidal zones, where they are limited to low-lying coastal areas with constant input of freshwater (Kathiresan & Bingham, 2001; Duke, 2017). Along the African Atlantic coast, the latitudinal extent of mangroves is constrained by two cold sea currents, the Benguela Current in the south, flowing northwards from South Africa to Angola, and the Canary Current flowing southwards from Morocco to Senegal (Saenger & Bellan, 1995) (Fig. 1). The presence and extent of mangrove stands in Africa are directly related to humidity, rainfall, and coastal relief and are known to have undergone extensive expansions and contractions in both range and abundance over time (Saenger & Bellan, 1995). Indeed, the study of mangrove cover is a primary source of information for modelling impacts of climate change (Alongi, 2015), as well as for palaeoclimatic reconstructions based on historical data from mangrove pollens (Bonefille et al., 1982; Leroy & Dupont, 1994; Scourse et al., 2005; Durogbo et al., 2010; Adeonipekun et al., 2016; Vallé et al., 2017).

Fig. 1
figure 1

Map of West and West-Central Africa with mangrove distribution (light green) and Aplocheilichthys spilauchen sampling sites (black dots) indicated. AC Angola Current, AR Atlantic Rise (red), BC Benguela Current, CC Canary Current, CVL Cameroon Volcanic Line (yellow bar), FP Freshwater plume (blue), GC Guinea Current, GGC Gulf of Guinea Current, SEC South Equatorial Current. Arrows indicate currents direction

Full size image

The Atlantic coast of Africa has an extremely narrow continental shelf, resulting in mangrove expansion during high sea level but contraction during low sea level (Saenger & Bellan, 1995; Scourse et al., 2005). In a high sea-level scenario, mangroves expand over new floodplains along the lower sections of rivers newly under tidal influence and salinity. In the opposite low sea-level scenario, fewer surfaces are exposed and mangroves recede to the newly reduced tidal zone, with the original stands replaced by other vegetation types. This is in contrast to regions with wide continental shelves, such as in Southeast Asia where significant mangrove expansion is related to low sea level, as shallow shelf areas previously covered by saltwater are exposed to tidal influence and freshwater inputs (Luther & Greenberg, 2009). A similar pattern is seen along the Brazilian coast, also with a wide continental shelf, where during periods of low sea level, many brackish water palaeochannels form and connect adjacent river drainages providing ideal conditions for mangrove expansion (Thomaz et al., 2015; Thomaz & Knowles, 2018). Similar to the described mangrove dynamics, the distribution of brackish water fishes is also highly dependent on freshwater inputs from effluent river systems and sea-level variations (Saenger & Bellan, 1995; Kathiresan & Bingham, 2001; Duke, 2017). Such alternating expansion and contraction dynamics can reveal interesting historical information when analysed in phylogenetic and temporal contexts, and this is particularly the case for taxa inhabiting regions known for major shifts in river outflows, climatic instability, sea level changes, and tectonic activity such as along the African Atlantic coast (Giresse, 2005; Goudie, 2005; Runge, 2007; Flugel et al., 2015).

Among the procatopodids, a cyprinodontoid group comprising mainly freshwater fishes endemic to Africa, the banded lampeye, Aplocheilichthys spilauchen (Duméril, 1861), is one of the few brackish water species and is a dominant species in mangrove habitats along the tropical and subtropical zones of the Atlantic coast (Fig. 1). Similar to other procatopodids and many other cyprinodontoids, A. spilauchen is known for a well-developed swimming capability (vagility) when compared to members of the Aplocheiloidei that in Africa are represented by the Nothobranchiidae (Huber, 1999). The latter is found mostly in shallow freshwater habitats, and some taxa even exhibit a seasonal life cycle (e.g. Nothobranchius), thus, presenting a completely different response to climatic and biogeographic processes. A river or an increase in the input of freshwater into an estuarine region likely represents range expansion opportunities for cyprinodontoid taxa, whereas for aplocheiloids, they more likely act as barriers (Bartáková et al., 2015).

Although A. spilauchen has also occasionally been collected from inland freshwater habitats, it primarily has a coastal distributional range, which closely mirrors the area covered by Atlantic coast mangrove forests, from the Senegal River to the Kwanza River in Angola (Wildekamp, 1995; Saenger & Bellan, 1995; Feka & Morrison, 2017) (Fig. 1). Aplocheilichthys spilauchen is currently recognized as the sole member of Aplocheilichthys and was the first African lampeye to be described (as Poecilia spilauchen), based on specimens from the Ogowe River in Gabon by Duméril (1861). It is readily distinguished from all other procatopodids by the presence of a combination of vertical grey banding along the flanks and dorsal, anal, and caudal fin in males and attainment of large body size (Fig. 2). Adults routinely reach up to 7.0 cm standard length while other procatopodids rarely exceed 4.0 cm standard length (Wildekamp, 1995). Likely, due to this distinctive pigmentation patterning and its presence in mangroves and other brackish water habitats, it has been assumed that A. spilauchen is a widely distributed species. Two other species, A. typus Bleeker, 1863 probably from Ghanaian coastal regions (Wildekamp, 1995), and A. bensonii (Peters, 1864) from Liberia were described but promptly synonymized with A. spilauchen by Günther (1866). A third species, A. tschiloangensis Ahl, 1928, described from the Tschiloango River in Cabinda, an Angolan territory north to the Congo River outlet, is now also considered a synonym of A. spilauchen (Huber, 1982; Wildekamp et al., 1986).

Fig. 2
figure 2

Males of A. spilauchen from localities along the Atlantic coast of Africa: a UFRJ 4150, 37.6 mm SL, Sangrougrou River, Senegal; b AMNH 275472, 51.0 mm SL, Farmoriah River, Guinea; c AMNH 59382, 52.8 mm SL, Freetown, Sierra Leone; d UFRJ 11487, 44.8 mm SL, Assimie, Ivory Coast; e MRAC 73–08-P-135–138, 41.0 mm SL, Ebrie Lagoon, Ivory Coast; f AMNH 226603, 46.4 mm SL, coastal Benin; g UFRJ 11484, 42.3 mm SL, coastal Nigeria (aquarium import); h MRAC 143378–387, 64.7 mm SL, Bioko Island, Equatorial Guinea; i MRAC 73–39-P-2168–179, 41.4 mm SL, mouth of Mirupururu River, Bioko Island, Equatorial Guinea; j AMNH 258331, 45.5 mm SL, mouth of Kouilou River, Republic of Congo; k AMNH 238526, 31.0 mm SL, Boma estuary, lower Congo River, Congo DRC; l SAIAB 84605, 40.0 mm SL, Dondo, Kwanza River, Angola

Full size image

Investigating seasonal variations in mangrove fish communities in Nigeria, Wright (1986) found that A. spilauchen was the dominant species, both in numbers and biomass and that its reproductive cycle was dependent upon freshwater inputs. Okyere (2012), in a study evaluating the use of A. spilauchen for mosquito larvae control along the Ghanaian coast, found that the species exhibited a low resilience to salinity increase, tolerating only a maximum of 4‰ salinity (Okyere, 2012). Both ecological studies suggest a dependence on freshwater in shaping the niche requirements and distribution of the species. Osteological information available for A. spilauchen reveals some inconsistencies regarding the configuration of the hypural plate supporting the caudal fin, usually an invariant feature among procatopodid species. The plate was considered as completely fused (Parenti, 1981), separated (Ghedotti, 2000), or bearing a small gap close to the compound caudal centrum (Costa, 2012).

Recently, in the first molecular dated analysis focused on the Procatopodidae, A. spilauchen was resolved as the sole member of a lineage sister to all other procatopodids, except Plataplochilus Ahl, 1928 (Bragança & Costa, 2019). Through the inclusion of two fossil calibrations in the sister families Valenciidae and Aphaniidae, it was estimated that A. spilauchen split from the remaining procatopodids in the early Miocene (around 23 mya). Considering the diversification patterns seen in all other procatopodid lineages, which intensified during the middle/late Miocene and Pliocene, the presumed lack of diversification within the A. spilauchen lineage is worthy of investigation.

Available ecological information (Wright, 1986; Okyere, 2012), conflicting osteological data (Parenti, 1981; Ghedotti, 1998, 2000; Costa, 2012) and a proposed early Miocene origin for the A. spilauchen lineage (Bragança & Costa, 2019) provide the impetus for the present study which aims to (1) investigate cryptic diversity within the A. spilauchen lineage through the application of both distance and coalescent species delimitation methods, (2) examine connectivity or genetic structuring between haplotypes, (3) estimate the divergence times for identified lineages, (4) perform ancestral area reconstructions, and (5) interpret the resulting temporal and biogeographic patterns in view of the current knowledge of the main historical events that have affected the African landscape and influenced Atlantic mangroves and coastal habitats.

Material and methods

Specimen preservation and fixation

Samples were collected using a variety of fishing techniques, and specimens were euthanized with clove oil or MS-222 anaesthetic solutions, in accordance with recommended guidelines for the use of fishes in research (Bennett et al., 2016). A small piece of muscle tissue or fin was taken from the right side of each specimen, or the entire fish was preserved in 95% ethanol in the field. The samples were stored at low temperatures at each of the following institutions: (1) NRF-South African Institute for Aquatic Biodiversity (NRF-SAIAB), Makhanda (Grahamstown), South Africa; (2) Ichthyological Collection of the Biology Institute of the Federal University of Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil; (3) Royal Museum for Central Africa (RMCA), Tervuren, Belgium; (4) American Museum of Natural History (AMNH), New York, USA; (5) National Museum of Natural History (USNM), Washington, USA and (6) Oregon State University (OS), Corvallis, USA.

Taxon sampling

Samples of newly generated and published sequences of the mitochondrial gene COI (Cytochrome c oxidase subunit 1) with a total of 679 bp from eight populations of A. spilauchen from across its geographical range were included in the analysis (Fig. 1). The present dataset includes samples from the (1) Kwanza and Lucala rivers in Angola (n = 7); (2) lower Congo River in the Democratic Republic of Congo (n = 2); (3) Kouilou River in the Republic of Congo (n = 2); (4) Ngounie River (Ogowe River drainage) and Komo river in Gabon (n = 3); (5) Ntem River in Equatorial Guinea (n = 3); (6) Aquarium import specimens from coastal Nigeria (n = 2); (7) Kakum River in Ghana (n = 3) and (8) Forecariah River at Conakry, Guinea (n = 2). Following Bragança & Costa (2019), three species were selected as outgroups: Plataplochilus miltotaenia Lambert, 1963 and P. pulcher Lambert, 1967, representing the first diverging lineage of the Procatopodidae and ‘Poropanchax’ normani (Ahl, 1928), representing the remaining procatopodid lineages sister to Aplocheilichthys. The generic name of ‘Poropanchax’ normani is between quotation mark to indicate that the species is not a member of Poropanchax sensu stricto as revealed in Bragança & Costa (2019), and awaits future taxonomic and nomenclatural resolution. A list of all included specimens and their respective catalogue numbers, localities, and GenBank Accession numbers are provided in Table 1. Other specimens, not included in the molecular analyses but photographed to illustrate pigmentation and body shape variation in different populations along the African coast, are shown in Fig. 2. The catalogue number and locality information of the specimens in Fig. 2 are listed in Online Resource 1.

Table 1 Species localities and Genbank and Bold Accession numbers: numbers in bold refer to sequences developed in the present study
Full size table

DNA extraction and sequencing

DNA was extracted from preserved tissues using the salting out method (Sunnucks et al., 1996), or the DNeasy Blood & Tissue Kit (Qiagen, Hilden) following the manufacturer’s standard protocol for animal tissue isolation. A fragment of the mitochondrial gene COI (Cytochrome c oxidase subunit 1) was amplified with universal primers LCO1490 and HCO 2198 and published protocols (Folmer et al., 1994). PCR products were purified with Exosap (Applied Biosystems), cycle sequenced using BigDye Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) and sequenced at the NRF-SAIAB using an ABI 3730xl DNA Analyzer (Applied Biosystems), or by Macrogen, South Korea.

Phylogenetic analyses, COI codon partitioning and evolution models

Phylogenetic analyses included newly generated and published sequences of the mitochondrial gene COI from A. spilauchen specimens from along the western coastal plains and mangroves of the Atlantic coast (see Table 1). Our COI alignment did not have frameshifts, indels or premature stop-codons, which are indicative of pseudogenes. The best models of evolution for each codon position for the maximum likelihood (ML) and Bayesian inference (BI) analyses were determined in PartitionFinder2 (Lanfear et al., 2016). The partitioned ML analysis was completed with GARLI (Zwickl, 2006), and bootstrap support was assessed from a majority rule consensus tree generated in Mesquite (Maddison & Maddison, 2019) from 1000 trees. The partitioned BI analysis was performed with MrBayes version 3.2 (Ronquist et al., 2012), and posterior probabilities were assessed with 5 million generations, sampling trees every 1000 generations. The first 25% of trees were discarded as burn-in. Both the ML and BI analyses were performed on the CIPRES Science Portal (Miller et al., 2010). The models of evolution implemented for each codon position in the ML and BI analyses were TRNEF + I + G, F81 + I, and TIM + G and SYM + G, F81 + I, and GTR + G, respectively.

For the species delimitation methods (described below), an ultrametric tree including only unique haplotypes (n = 15) as required by the GMYC method was performed in Beast 1.8.2 (Drummond et al., 2012), and the parameters were defined as follows: an uncorrelated relaxed clock model with a lognormal distribution (Drummond, et al., 2006), with 10 million generations with a sampling frequency of 1000. The value of parameters of the analyses, convergence of the MCMC chains, sample size and the stationary phase of chains were evaluated using Tracer 1.6 (Rambaut et al., 2014). A Birth and Death Incomplete Sampling speciation process for the tree prior (Stadler, 2009), indicated for datasets with incomplete sampling, was used, and the analysis included only ‘Poropanchax’ normani as outgroup. For this analysis, the evolution model HKY + G was inferred in JModeltest (Darriba et al., 2012) for the entire COI gene. Saturation levels were checked in Dambe5 (Xia, 2013), according to the algorithm proposed by Xia et al. (2003).

Species delimitation

Four species delimitation methods were applied to investigate the diversity within A. spilauchen. The traditional barcoding genetic distance method (GD) (Herbert et al., 2003) and the Automatic Barcode Gap Discovery (ABGD) (Puillandre et al., 2012), both relying on genetic distances between haplotypes, and the General Mixed Yule Coalescent (GMYC) (Fujisawa & Barraclough, 2013) and Bayesian implementation of the Poisson Tree Processes (bPTP) (Zhang et al., 2013), both coalescence approaches.

For the GD method, we used the Kimura-2-parameters model (K2P) (Kimura, 1980) to estimate the pairwise genetic distances between the different A. spilauchen populations in MEGA 7 software (Kumar et al., 2016). Haplotypes with a genetic divergence higher than 3% were considered as belonging to distinct operational taxonomic units (OTUs) following the expected pattern for genetic divergence between fish species (Ward, 2009). ABGD also relies on genetic distances to identify the threshold between interspecific (speciation) and intraspecific (populational) processes within the dataset (Puillandre et al., 2012). The main difference between GD and ABGD is that the ABGD allows a more refined search; once the estimated genetic divergence between groups (putative species) is calculated, it is recursively applied to the previously delimited groups until no more putative OTUs are recognized. The ABGD result is presented relative to a spectrum of P values (prior intraspecific values) in which a 0.001 value assumes a minimum intraspecific variability, and a 0.1 value assumes a maximum intraspecific variability. The different results relative to each significant P value intervals are presented, and to identify which P value interval best describes the diversity within the dataset, a congruence with other methods is expected and/or a support from traditional morphological alpha taxonomy. The ABGD analysis was performed in the ABGD server website (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html) following the default parameters.

The coalescence species delimitation method, GMYC (Fujisawa & Barraclough, 2013), requires an ultrametric tree to distinguish between intraspecific (coalescent) and interspecific (speciation) threshold patterns, because it relies on branch length differences to define OTUs. The bPTP coalescent method (Zhang et al., 2013), in contrast to GMYC, relies on the number of nucleotide substitutions between haplotypes to resolve intraspecific and interspecific patterns. When performing bPTP, two results are provided, the BI and the ML solutions, and in each, the recognized OTUs are depicted graphically in a tree. The same reduced dataset tree as for GMYC was used when performing bPTP. The GMYC analysis was performed on the Exelixis Lab's server https://species.h-its.org/gmyc/ following default parameters, and the bPTP was performed on the Exelixis Lab's server https://species.h-its.org/ptp/ following default parameters except for a 20% burn-in.

Haplotype network

The haplotype network was inferred, using the Minimum Spanning Network method (Kruskal, 1956) in PopART 1.7 (Leigh & Bryant, 2015), including only A. spilauchen haplotypes. The following eight geographic areas/drainages were delimited: (A = Gabon, B = Angola, C = Guinea, D = Kouilou, E = Lower Congo, F = Nigeria, G = Ghana, H = Equatorial Guinea). These areas were delimited based on the presence or absence of connectivity between haplotypes as suggested by the genetic divergence and species delimitation results. Therefore, the area ‘Angola’ includes all haplotypes from the Kwanza and Lucala rivers, the latter being one of the main tributaries of the lower Kwanza River. The area ‘Gabon’ includes all haplotypes from the Ngounie River, a tributary of the lower Ogowe River, and from the Komo River.

Time-calibrated analysis

The time-calibrated analysis was constructed in BEAST v.1.8.2 (Drummond et al., 2012), including all sequences available, and an uncorrelated relaxed clock model with a lognormal distribution (Drummond et al., 2006). Bayesian Inference was performed with 50 million generations with a sampling frequency of 1000. The values of parameters of the analysis, convergence of the MCMC chains, sample size, and the stationary phase of chains were evaluated using Tracer 1.6 (Rambaut et al., 2014). The same Birth and Death Incomplete Sampling speciation process was selected for the tree prior (Stadler, 2009). Since no fossil procatopodids are currently known three relevant calibration points based on the age estimates from Bragança & Costa (2019) were selected. These secondary calibrations correspond to (1) the Procatopodidae basal node (prior setting: normal distribution, mean = 30.1, standard deviation = 1.3); (2) the node representing the split between the A. spilauchen clade and ‘Poropanchax’ normani (prior setting: normal distribution, mean = 23.1, standard deviation = 1.5); and (3) the node representing Plataplochilus diversification (prior setting: normal distribution, mean = 3.2, standard deviation = 0.9). A normal distribution prior was selected for the secondary calibration points as recommended by Ho (2007) and Ho & Phillips (2009). At the end of the analysis, a burn-in of 25% of the retained trees was performed in TreeAnnotator.

Ancestral area reconstruction

The ancestral area reconstruction method, S-DIVA (Statistical Dispersal-Vicariance Analysis) (Yu et al., 2010), was performed in RASP 3.2 (Yu et al., 2015) to reconstruct the expected ancestral areas during A. spilauchen diversification. S-DIVA is an event-based method that incorporates statistical uncertainty into both phylogenetic and ancestral state reconstructions and assumes vicariance as the most probable event to explain biogeographical evolution within a group, attributing higher costs for dispersal and extinction events. For the S-DIVA analysis, both extinctions and reconstructions were considered, and at each node, the frequency of alternative ancestral area reconstructions generated for each tree in the dataset is shown. The analysis was performed on the resulting trees from the time-calibrated analysis from Beast 1.8.2 (Drummond et al., 2012), and the condensed tree was defined as the summary tree from that same analysis. Prior to the analysis, the outgroups were removed in RASP 3.2, and a post-burning of 12,500 trees was carried out. The maximum number of areas included in ancestral distributions was set to 3, but based on current knowledge from ecology, distribution, and genetic divergence between populations, only ancestral distributions between contiguous areas were allowed.

Results

Phylogenetic analysis (Online Resource 2)

ML and the BI analyses recovered the same topology with most nodes receiving posterior probability values > 95 and bootstrap values > 80, except for nodes representing the most recent divergences (Online Resource 2). An early divergence of Plataplochilus species and a sister group relationship between the A. spilauchen lineage and ‘Poropanchax’ normani were recovered with strong support. Within Aplocheilichthys, a sister group relationship between haplotypes from Ghana and Guinea was supported with a high posterior probability but a comparatively lower bootstrap value (Online Resource 2). Haplotypes from Gabon were resolved as sister to all remaining haplotypes with high support. Another node, also strongly supported, represents the split between Nigerian + Equatorial Guinean haplotypes from those from Kouilou, lower Congo, and Angola. A sister group relationship between Nigerian and Equatorial Guinean haplotypes was supported by maximum values in both analyses. The node-grouping haplotypes from Kouilou, lower Congo and Angola were supported with high posterior probability but a comparatively lower bootstrap value. However, within that clade, all nodes received lower support values, except the node representing the split between Kouilou haplotypes and those from lower Congo and Angola, which received strong support.

Species delimitation (Fig. 3)

Fig. 3
figure 3

Ultrametric tree and OTUs delimited by each species delimitation method (represented by different colours). Numbers indicate posterior probability values

Full size image

Both genetic distance methods provided similar results, but as expected, with a decrease in intraspecific variability value (P value) in ABGD more OTUs were delineated. The ABDG result with a P value between 0.0129 and 0.0599 recovered the same six OTUs as delimited by the GD method: (1) Ghana, (2) Guinea, (3) Gabon, (4) Nigeria, (5) Equatorial Guinea, (6) Kouilou + lower Congo + Angola, whereas one more OTU was identified between P values of 0.0046 and 0.0077, with the Kouilou lineage considered as distinct from a lower Congo + Angola OTU. The absolute pairwise genetic distances for the OTUs identified based on the GD method range from 8 to 22%, which are considerably higher than the heuristic standard threshold of 2–3% commonly applied to recognize species limits (Table 2). The highest genetic distance (22%) was between specimens from Guinea and Equatorial Guinea while the lowest (2%) was recorded for specimens from Kouilou, lower Congo and Angola.

Table 2 Genetic distance within A. spilauchen haplotypes conducted using the K2-parameter model in MEGA 7
Full size table

The Bayesian ultrametric tree (Fig. 3) had the same topology as the phylogenetic analyses (Online Resource 2), with most nodes supporting putative OTUs with maximum posterior probability values. There was concordance of the six OTUs defined by ABGD (P value between 0.0129 and 0.0599), GD, GMYC and the bPTP ML solution result, while bPTP BI solution returned the same seven OTUs as ABGD when assuming the lowest P values. We considered the six OTUs defined by the traditional barcoding, the ABGD including P values between 0.0129 and 0.0599, the GMYC and the bPTP ML solution as putative species but defer making any taxonomic changes pending a detailed osteological and morphometric study (see Discussion).

Haplotype network (Fig. 4b)

Fig. 4
figure 4

a Map showing the sampled localities and the male colouration in life at each site; b mitochondrial gene COI Minimum Spanning haplotype Network for A. spilauchen haplotypes. The number of nucleotide changes is indicated in parentheses

Full size image

A minimum spanning network portraying genealogical relationships among haplotypes revealed considerable divergence between geographically separated populations within A. spilauchen. In most cases, the number of nucleotide substitutions separating populations was extremely high ranging from a minimum of 43 up to 133. The only exception to this being among haplotypes from the neighbouring Kouilou, lower Congo, and Angolan populations which form a cluster with a maximum of 12 nucleotide changes separating Kouilou haplotypes from the Angolan and lower Congo haplotypes.

Time-calibrated analysis (Fig. 5)

Fig. 5
figure 5

Time-calibrated phylogeny of the Procatopodidae and A. spilauchen obtained from the Bayesian dating analysis in BEASTv.1.8. Bars represent maximum and minimum date estimates for each node (values in brackets), and the numbers are nodes divergence mean ages. The colours in the time bar are a reference to the proposed time extension of main palaeogeographic and palaeoclimatic events during A. spilauchen evolution: the green bar represents the Miocene Climatic Optimum; the red bar represents the Miocene Climatic Transition; the orange bar represents the late Miocene aridification; the brown bar represents the Pliocene–Pleistocene climatic instability; the blue bar represents the onset of the modern Congo River outlet and the black bar represents an increase in the volcanic activity of the CVL

Full size image

As expected, given a single marker dataset and uncertainties inherent in the use of secondary calibrations, most nodes on our time tree have 95% HPD values spanning large time intervals, and there is no question that considerable uncertainty persists regarding the precision of the resultant dating scheme. Nonetheless, despite these caveats, our scheme does provide a broad estimate of comparative divergence times among geographically disparate populations of A. spilauchen (Fig. 5). Based on the current analysis, the onset of diversification within the A. spilauchen lineage began during the Middle Miocene (14.8 mya, 95% HPD 8.6–20.7 mya), represented by the split of a west African lineage comprising specimens from Ghana and Guinea. These two lineages subsequently diverged during the late Pliocene (3.5 mya, 95% HPD 0.6–8.7 mya). A split between specimens from Gabon and all the remaining populations occurred during the late Miocene (9.5 mya, 95% HPD 4.8–14.9 mya). This was followed by another in the transition between the late Miocene and the Pliocene (5.6 mya, 95% HPD 2.3–9.5 mya), in which two clades diverged, one including specimens from Nigeria and Equatorial Guinea and the other those from the Kouilou, lower Congo and Angola. Later, the lineage including Nigerian and Equatorial Guinean specimens diverged at the onset of the Pleistocene (2.5 mya, 95% HPD 0.8–5.0 mya). The remaining lineage diversification time estimates are extremely recent ranging from around 1 mya to the near present.

Ancestral area reconstruction (Fig. 6)

Fig. 6
figure 6

S-DIVA Ancestral area reconstructions of A. spilauchen, and input chronogram resultant from BEASTv.1.8 dated analysis. Colours for each designated area are presented in legend box

Full size image

Biogeographic ancestral area reconstruction inferred the current distribution pattern of A. spilauchen lineages to have likely resulted from repeated dispersal and vicariance events (Fig. 6). S-DIVA postulates 6 dispersal, 6 vicariance, and 1 extinction event. The root node, corresponding to the onset of diversification of the A. spilauchen lineage, estimated the areas FGH (Ghana, Nigeria and Equatorial Guinea) and CFG (Guinea, Ghana, and Nigeria), as equally most probable ancestral areas (49% probability). From this root node, a vicariance event is suggested between areas AFH (Gabon, Nigeria, and Equatorial Guinea) and CG (Guinea and Ghana). At the node corresponding to the split between Guinean and Ghanaian haplotypes, the ancestral area CG was recovered (96% probability), and a vicariance event was suggested as the cause of the split between the two areas. The ancestral area AFH was identified for the node in which the Gabon haplotypes split from the remaining haplotypes (95% probability), and it is possible to identify a duplication event (within area speciation) in area A (Gabon) followed by a vicariance event in which one lineage is restricted to Gabon and the other present in area AFH. The area AFH was considered as the ancestral area for the node in which Nigeria and Equatorial Guinea haplotypes split from the Kouilou, lower Congo and Angolan haplotypes (93% probability). Vicariance was recovered as the explanation for the split between areas FH (Nigeria and Equatorial Guinea) and ADE (Gabon, Kouilou, and lower Congo). The ancestral area FH was estimated for the node in which Nigeria and Equatorial Guinea haplotypes split (96% probability), caused by a vicariance event. For the node which corresponds to the split between Kouilou haplotypes and Angola + lower Congo haplotypes, the ancestral area ADE (Gabon, Kouilou, and lower Congo) (93% probability) was estimated, followed by an extinction event in area A (Gabon), a dispersion to area B (Angola) and finally a vicariance event between areas BE (Angola and lower Congo) and D (Kouilou). For the split between Angola and lower Congo haplotypes, the ancestral area BE was delimited (100% probability) and considered to be the result of a vicariance event.

Discussion

Prior to the current study, A. spilauchen had been considered to be widely distributed with a range extending along much of coastal west and west-central Africa. This was based on the assumption that a higher salinity tolerance, relative to that of other procatopodids, would have allowed the species to maintain population connectivity across this extensive geographical range. In addition, an apparent lack of variability in pigmentation patterning between individuals from geographically disparate populations supported a widespread species hypothesis, which in turn resulted in no further investigation of morphological differences between them. Unfortunately, the ongoing COVID-19 pandemic has prevented us from undertaking a morphological investigation to accompany the present study. Loan of museum materials is currently not possible, and we are unable to examine the type specimens of three previously synonymized taxa, A. typus (Ghana), A. bensonii (Liberia) and A. tschiloangensis (Cabinda, Angola), or investigate potential osteological (Parenti, 1981; Ghedotti, 2000; Costa, 2012) or morphometric variation among populations. Further investigation of potential phenotypic differentiation between the molecular lineages identified here is necessary prior to formalizing any taxonomic conclusions, and consequently, we must defer such actions to a future contribution.

The present study has, however, uncovered considerable structuring and genetic divergences between populations that are frequently 4–11 times higher than the traditionally employed sequence divergence heuristic threshold of 2–3% for teleostean conspecifics (e.g. Pereira et al., 2013; Decru et al., 2016; Iyiola et al., 2018; Arroyave et al., 2019) (Fig. 4b, Table 2). Here, we focus our discussion on the potential drivers of diversification and mechanisms that are likely to have shaped the contemporary distributions of lineages within this complex.

The onset of A. spilauchen diversification coincides with one of the main climatic shifts during the Neogene (Fig. 5). Initially, a period known as the Middle Miocene Optimum (between 16 and 14.8 mya) was warm and humid with high precipitation. This was followed by the Middle Miocene Climatic Transition (14.8–12.9 mya) which was much dryer leading to a significant drop in sea levels and increased aridification across the African continent (Flower & Kennett, 1994; Kender et al., 2009, 2014; Herold et al., 2011; Frigola et al., 2018). During the Middle Miocene Optimum, a moist and warm climate promoted the expansion of tropical vegetation, which covered most of the continent expanding even to high latitudes (Lovett, 1993; Maley, 1996; Plana, 2004; Stanley et al., 2005; Frigola et al., 2018). During this time, sea level rise likely engendered mangrove expansion throughout the Atlantic coast, and consequently potential expansion of an ancestral range for A. spilauchen. Our ancestral area analysis estimates an expansion from an area that corresponds to the present-day Ghanaian, Nigerian and Equatorial Guinean coastlines to adjacent Gabon to the south and Guinea to the west (Fig. 6). During the Middle Miocene Climatic Transition, this area experienced a marked contraction of mangrove habitats isolating Aplocheilichthys lineages in the westernmost coastal regions of Guinea and Ghana, and another in Nigeria, Equatorial Guinea and Gabon (Fig. 6).

Studies on the evolution of the Congo River drainage (Beadle, 1981; Burke, 1996; Giresse, 2005; Goudie, 2005; Stankiewicz & de Wit, 2006; Runge, 2007), deposition patterns in the Congo deep sea fan (Lavier et al, 2001; Leturmy et al., 2003; Lucazeau et al., 2003; Anka & Séranne, 2004; Anka et al., 2009; Savoye et al., 2009) and freshwater fish diversification in the basin (Goodier et al., 2011; Schwarzer et al., 2011; Alter et al., 2015, 2017; Arroyave et al., 2020; Stiassny & Alter, in press) reveal a complex geologic history for the most diverse river basin in Africa. In the western basin, sedimentary studies suggest a protracted history of shifting and intermittent outflow of the Congo River into the Atlantic during the Cenozoic. Although a final consensus has yet to be reached, a single high-energy capture event is now generally considered to have established the current Congo outlet to the Atlantic shortly after the Miocene-Pliocene transition. A Pliocene capture (5.3–2.6 mya) is supported by an increase in sediment deposition in the Congo fan and by tectonic activity along the Atlantic Rise during that time (Lavier et al., 2001; Leturmy et al., 2003; Lucazeau et al., 2003; Anka & Séranne, 2004; Anka et al., 2009; Savoye et al., 2009). High humidity during this period is also suggested by palynological data and probably these events facilitated a range expansion from Nigeria, Equatorial Guinea, and Gabon southwards to the regions in and around the newly formed lower Congo estuary at the border between the Democratic Republic of Congo and Angola (Fig. 6). Recent data on the extent of the Congo River plume and its influence on sea surface salinity and temperature over expansive coastal regions (Materia et al., 2012; Denamiel et al., 2013; Chao et al., 2015) highlights the importance of the origin of the present-day Congo River outlet for the expansion of mangroves throughout the region. During wet seasons, the Congo plume connects with the Niger River plume, and a lower salinity is recorded for the entire Gulf of Guinea up to and including the region around the Congo outlet.

West African aridification during the Pliocene (5.3–2.6 mya) is documented by both palynological and sediment data sampled from ocean drilling sites along the coast (Bonnefille et al., 1982; Leroy & Dupont, 1994; Vallé et al., 2017). Study of sedimentological sequences reveals a significant reduction in river discharge and an increase in wind borne grass pollen, both indicative of an arid climate. In addition, a concomitant reduction in mangrove pollen (Rhizophora spp.) suggests that aridification resulted in mangrove reduction and fragmentation. The late Pliocene aridification likely facilitated two vicariance events observed in our study: one between the Ghanaian and Guinean lineages and the other between the Nigerian and Equatorial Guinean lineages. Studies at sites close to Cape Blanc in Mauritania and just south of the mouth of the Senegal River found similar patterns, indicating an onset of aridification around 3.4 mya (Bonnefille et al., 1982; Leroy & Dupont, 1994). A similar study at a site near the Comoé River outlet, identified aridification in north western Africa between 3.5 and 2.9 mya (Vallé et al., 2017), a timing consistent with the proposed date for the split between Ghanaian and Guinean lineages (3.5 mya) estimated in the present study (Fig. 5)).

Other sites with detailed palynological information are near the Niger Delta, where a pronounced aridification with an increase in wind borne grass pollen and marked reduction of Rhizophora pollen is recorded between 2.7 and 2.0 mya (Durugbo et al., 2010; Adeonipekun et al., 2016). These dates are also consistent with the findings of the present study where the split between Nigerian and Equatorial Guinean lineages is estimated to have occurred around 2.5 mya (95% HPD 0.8–5.0 mya) (Fig. 5). Vicariance between Nigerian and Equatorial Guinean lineages may also be related to increased activity of the Cameroon Volcanic Line, a chain of mountains and volcanos that extends from the Cameroonian Highlands to the volcanic islands in the Gulf of Guinea (Deruelle et al., 1991; Marzoli et al., 2000; Burke, 2001). Bioko Island is located about 60 km from the mainland and experienced volcanic activity during the same period with volcanism of Cameroon and Manengouba mounts (around 3.0–1.0 mya) near the coast (Deruelle et al., 1991; Marzoli et al., 2000). Volcanism during low sea level would have had major impacts on water chemistry, temperature and coastal sedimentation directly impacting mangrove cover throughout the region. Recently, a comprehensive study of population connectivity of the mangrove, Rhizophora racemosa, around the Cameroonian and Equatorial Guinean coastline found high levels of genetic structuring related to the volcanic activity of Bioko Island and the Bioko-Cameroon land bridge formation (Ngeve et al., 2016). Here, we suggest that volcanism and the land bridge formation likely also affected A. spilauchen and resulted in the disruption of connectivity (gene flow) between Nigerian and Equatorial Guinean lineages.

The late Pleistocene-Holocene is known for considerable climate instability (Marius & Lucas, 1991; Scourse et al., 2005; Malounguila-Nganga et al., 2017; Maley et al., 2017; Molliex et al., 2019). Although our sampling of A. spilauchen from the southern extent of its range is limited, ancestral area reconstruction posits two additional vicariance events among haplotypes in the Kouilou, lower Congo and Angolan cluster during this timeframe (Fig. 6). Our results indicate a recent divergence and differentiation of the Kouilou haplotypes from the lower Congo and Angolan cluster. However, the low levels of genetic divergence detected, when compared to the significant genetic divergence between the other A. spilauchen populations, suggest the possibility of a continued gene flow between all of these populations (Fig. 4b, Table 2). A recent, and ongoing, connection is suggested by the influence of the Congo River freshwater plume, which is periodically carried south by the Angola current (Kopte et al., 2017), reaching beyond the Nyanga and Kwanza river ouflows (Denamiel et al., 2013). Similarly, an increase in mangrove vegetation persisting since the last deglaciation is evidenced by pollen data from the mouth of Congo (Scourse et al., 2005). Considerably denser population sampling of A. spilauchen throughout this southern region will be necessary to resolve this issue.

Conclusion

Despite relatively limited sampling, the applications of species delimitation, phylogeographic, and phylogenetic methods reveal a pattern of genetic differentiation within A. spilauchen that is concordant with a series of historical events recorded since the Middle Miocene. Our study highlights extremely low connectivity between most populations and a time-calibrated phylogeographic pattern that lends support to the novel hypothesis that a major driver of diversification within the lineage has been the shifting dynamics of coastal mangrove forest cover over time. We report, for the first time, a pattern of diversification within a lineage of brackish water fish that is concordant with the historical distribution of coastal mangroves forests, the predominant brackish water habitat of the focal species throughout its range.

From a conservation perspective, these results are of considerable significance since many brackish environments, particularly the mangrove forests, are highly threatened by coastal development and the exploration and extraction of hydrocarbons along the African Atlantic coastline (Alongi, 2015; Feka & Morrison, 2017). It has been estimated that mangroves are disappearing at a rate of 1–2% per year, suggesting the complete disappearance of these societally and biologically important ecosystems within the century (Alongi, 2015). The unexpected diversity observed in the A. spilauchen complex suggests that many of the other taxa that share a similar distribution associated with the same dynamic mangrove habitats may hide a significant, but currently undocumented, diversity under threat. Before formal taxonomic and nomenclatural changes can be undertaken, morphological analyses of the distinct populations (OTUs) identified herein are needed, and conservation assessments for each are critical as many of these potential new species are likely highly threatened by ongoing coastal development throughout the region.