Genetic and morphological variation of Woodland Kingfisher Halcyon senegalensis reveals cryptic mitochondrial lineages and patterns of mitochondrial–nuclear discordance

The Woodland Kingfisher Halcyon senegalensis is widely distributed throughout sub-Saharan Africa and occupies a wide variety of woodland and savannah habitat. Thus far, three subspecies have been described based on morphological variation. In the present study, using western, eastern and southern African populations, we examined the relationship between morphological and genetic divergence among two named subspecies, H. s. cyanoleuca and H. s. senegalensis, using three mitochondrial markers (CO1, Cytb, 16S) and two nuclear markers (FIB5 and RAG1). Southern birds showed clear evidence for morphological divergence, with a longer wing and tail length, when compared with eastern and western birds. Phylogenetic analyses using Bayesian methods identified two well-characterised genetic clusters, representing the two subspecies. We determined that H. s. senegalensis and H. s. cyanoleuca are closely related subspecies that split recently, approximately 0.66–1.31 MYA in the Pleistocene. Furthermore, genetic substructure was evident within H. s. senegalensis, with three distinct genetic clusters in each region. The separation between the Ghana+Gabon and Uganda lineages of H. s. senegalensis occurred approximately 0.12–0.57 MYA. Nuclear–mitochondrial discordance was detected, however, wherein the pattern of divergence was not detected in the RAG1 and FIB5 sequences. Our results suggest that climate change, biogeographic barriers and local adaptation has played a role in the diversification of Woodland Kingfishers in Africa.

as geographic isolation by distance or barriers (Wright 1943;Dobzhansky and Dobzhansky 1970;Greenwood 1980;Matthiopoulos et al. 2005;Taylor and Friesen 2012). The level of differentiation is dependent on the balance of gene flow, genetic drift and natural selection along environmental gradients (Rice and Hostert 1993;Avise 2004;Pinho and Hey 2010). The scale of differentiation varies greatly among species and may result in speciation (Slatkin 1987). In addition, widely distributed species may be further subdivided into subspecies. Species display phenotypic distinctiveness through to reproductive incompatibility (de Quieroz 2007), whereas subspecies are reproductively compatible but display at least one heritable trait and may be associated with a specific geographic area or ecological niche (Wallin et al. 2017). Accurately determining processes involved in speciation and differentiation is challenging and complex and has not been fully explored for several avian species on the African continent. Kingfishers are an example of a widely distributed species that is understudied.
Kingfishers (Alcedinidae) consist of three subfamilies (Halcyoninae or Daceloninae, Alcedininae, and Cerylinae), 17 genera and 92 species, of which 18 species occur in Africa (Fry et al. 1988). Thus far, studies on the phylogenetic analysis of kingfishers are limited, with disagreement regarding the basal family based on several lines of evidence (Maurer and Raikow 1981;Fry et al. 1988;Sibley and Ahlquist 1990;Johansson and Ericson 2003;Moyle 2006). The Woodland Kingfisher Halcyon senegalensis (Linnaeus, 1766) belongs to the subfamily Halcyoninae. However, uncertainty remains regarding the relationships within the subfamily, especially at the base of the radiation (Moyle 2006). The Woodland Kingfisher is widely distributed in sub-Saharan Africa (Moyle 2006), occupying a wide variety of woodland habitats. Although the species faces numerous threats, the International Union for Conservation of Nature has assessed it as Least Concern because of its extremely large distributional range (~20 100 km 2 ) and stable population trend (Birdlife International 2021). Three subspecies are currently recognised (Figure 1): Halcyon senegalensis fuscopileus Reichenow, 1906 occurs from Sierra Leone to southern Nigeria and south to the Democratic Republic of Congo (DRC) and northern Angola; Halcyon senegalensis cyanoleuca (Vieillot, 1818) is found from southern Angola and western Tanzania to South Africa; and Halcyon senegalensis senegalensis (Linnaeus, 1766) is distributed in Senegal and Gambia to Ethiopia and northern Tanzania. It has been reported that H. s. fuscopileus populations are found in forest habitat and are mainly or entirely resident, whereas H. s. cyanoleuca and H. s. senegalensis occur in well-developed woodland, such as in tall Acacia stands and Mopane, and are largely migratory (Fry et al. 1988). All three subspecies differ morphologically: H. s. fuscopileus is smaller and has a crown that is dark brownish grey, with a mantle and breast greyer than other subspecies; H. s. senegalensis and H. s. cyanoleuca are similar in colour, with H. s. cyanoleuca having a dark strip running behind the eyes (Fry et al. 1988). The migratory movement patterns, drivers of migration, phenology and phylogeography of intra-African migrant species like the Woodland Kingfisher are only just being understood. Available knowledge indicates trans-equatorial migration in H. s. cyanoleuca, with breeding grounds in southern latitudes of South Africa (reportedly >23°S), although resident populations are reported in Tanzania; H. s. senegalensis has resident populations on both sides of the Equator, with migrants of this subspecies making northern non-breeding movements into sub-Saharan Africa; migrant H. s. cyanoleuca has been recorded in the breeding ranges of H. s. senegalensis, however no breeding is reported to overlap, and the species is described as monogamous.
To date, limited phylogenetic studies have been completed for Woodland Kingfisher (Maurer and Raikow 1981;Fry et al. 1988;Sibley and Ahlquist 1990;Johansson and Ericson 2003;Moyle 2006) and no genetic studies have been conducted to elucidate the phylogenetic relationships of H. s. cyanoleuca and H. s. senegalensis. Modern molecular tools can be used to provide insight into the evolutionary history of taxa (Jetz et al. 2012) and can be used for the discovery of previously unrecognised cryptic species (Bickford et al. 2007). This study examined two of the Woodland Kingfisher subspecies using genetic methods, namely sequencing of mitochondrial DNA (mtDNA) and nuclear genes. mtDNA has been preferentially used for taxonomic and phylogeographic studies because of its lack of recombination, ease of amplification, and maternal inheritance (Boonseub et al. 2009). Sequence-based nuclear DNA variation can be used to complement mtDNA markers to further understand systematic relationships and genetic variation within a species (Lessa 1992). Thus, a combination of mtDNA and nuclear markers is recommended in demographic studies, as there are several factors (e.g. non-neutrality, extreme rate variation and recombination) that could confound inferences provided by mtDNA alone.
Periodic climatic oscillations between wet and dry conditions have occurred over the last 4-5 million years (Zachos et al. 2001) and have led to significant habitat modification or fragmentation, resulting in diversification of species (Hewitt 2003). However, the impact of climatic fluctuations varies per taxa as well as per geographic range (Stewart et al. 2010). Therefore, we investigated the phylogenetic relationships and genetic history of divergence, for H. s. senegalensis lineages from West Africa (Ghana) and East Africa (Uganda), and H. s. cyanoleuca from South Africa, using mitochondrial and nuclear molecular markers. We hypothesise that the Woodland Kingfisher in Africa has been subjected to differentiation or speciation owing to historic isolation and subsequent connectivity between populations as suitable habitat for the species expanded and contracted across Africa.

Sample collection
Across western, eastern and southern African, between November 2015 and January 2019, Woodland Kingfishers were trapped using varying numbers and lengths of mist nets, as well as spring traps baited with superworms Zophobas morio. The traps were deployed during morning (06h00-10h00) and evening (15h00-18h00) sessions. Field identification of subspecies during sample collection was based on a subspecies known range as well as the diagnostic eye stripe of H. s. cyanoleuca and the darker plumage of H. s. fuscopileus (H. s. fuscopileus was not encountered during the study). A total of 60 individual Woodland Kingfishers were trapped and sampled at least once, with 15 additional re-traps that occurred in southern Africa.
In Ghana (western Africa), blood samples were collected from ten H. s. senegalensis individuals in the Greater Accra Province (n = 5) and in the Central Province (n = 5), which are ~150 km apart. In Uganda ( Woodall 2001, cited in del Hoyo et al. 2014 individuals trapped in the Central Region. In South Africa (southern Africa), blood samples were collected from 38 H. s. cyanoleuca individuals trapped in Limpopo Province. Estimated distances between the focal countries are as follows: Ghana and Uganda, ~5 000 km; Ghana and South Africa, ~7 000 km; and Uganda and South Africa, ~3 700 km. Generally, samples were collected during the breeding season across the subregions. Samples were collected in Ghana between June and September, in Uganda between July and August, and in South Africa between November and January. Blood samples were collected by a brachial wing venipuncture, using a 27-gauge needle and 100 µl capillary tubes, and were stored in lysis buffer (Seutin et al. 1991). All biological materials collected were stored at the Biobank of the South African National Biodiversity Institute (SANBI).
Trapped birds were ringed using individually coded aluminium rings which followed the ringing scheme in each country, as well as with plastic colour rings that followed a unique combination. The metal rings ensured individuality of samples as well as the identification of individual birds if re-trapped, whereas the coloured plastic rings facilitated identification of the individuals when free and perched. All trapped individuals were weighed and measured before the blood tissue samples were collected, and then were released immediately after sampling.

Morphometric measurements
The morphometric measurements made on trapped H. s. cyanoleuca and H. s. fuscopileus were as follows: • Mass (g) -in the first two years of the study this was measured using a spring balance (Pesola 20100 Micro-Line Metric Spring Scale, 100g), whereby the bird is weighed in a bird bag and the bird bag alone was then weighed to determine the mass of the bird; in subsequent years, a digital scale (Pesola PPS200 Professional Digital Pocket Scale, 200g) was used, during which a small plastic container was tared before each measurement and the bird was placed in the container and then weighed. • Wing length (mm) -measured using a wing rule, placing the bend of the wing against the top of the rule, flattening the wings and feathers so that the measure is maximised, and taking the reading from the tip of the longest wing feather (the primaries). • Tail length (mm) -measured using a flat rule and taken from the base of the tail to the tip of the longest tail feather. • Head length (mm) -measured using a digital calliper (0-150 mm) and taken from the back of the skull to the front of the skull; this measure excluded the length of the culmen from the total head length. • Tarsus length (mm) -measured using a digital calliper (0-150 mm) and taken from the notch on the metatarsus (where it meets the tibiotarsus) to the top of the bone above the bent toes.

Multivariate analysis
Summary statistics of variation among the morphometric variables were estimated with PAST software (Hammer and Harper 2005). Principal components analysis (PCA) was performed to assess significant variation of five morphometric characters among the localities and between the subspecies, for the different sexes, also using PAST software (Hammer and Harper 2005). Four of the morphometric characters were direct measurements for each bird, namely: head length, tarsus length, wing length and tail length, whereas the body mass index (BMI) was calculated as a relationship of tarsus length to mass. This was done using the following formula: BMI = mass (g) / tarsus length (mm) (Nesbitt et al. 2008). All morphometric measurements were taken using the same devices to reduce the effects of measurement errors. Furthermore, images of the birds and the sampling process were taken for verification of field data and notes. Measurement error was tested using one-way analysis of variance (ANOVA) with Kruskal-Wallis tests (Kruskal and Wallis 1952), by comparing measurements taken by the samplers for each variable. The Kruskal-Wallis test compares variation between two or more independent samples of equal or different sample sizes. To minimise size-related dissimilarity, a PCA with log-transformed morphometric variables was also performed for data on adult birds, using PAST software. Log-transformation of continuous datasets removes the skewness to make data as normal as possible, to ensure more valid analysis. Multivariate analysis of variance (MANOVA) was also conducted in PAST to test whether there was significant morphometric variation between sexes, for the three localities and two subspecies.
All reactions were conducted using Ampliqon Taq Master Mix RED in a final reaction volume of 15 µl, containing 6.25 µl of Ampliqon Taq Master Mix RED (Ampliqon A/S, Odense, Denmark) (0.1 M Tris/HCl, pH 8.5, (NH 4 ) 2 SO 4 , 4 mM MgCl 2 , 0.2% Tween 20, 0.4 mM deoxynucleotides, 0.2 units µl -1 Taq DNA Polymerase, inert red dye and stabiliser), 0.1 µM of each the forward and reverse primers (Thermo Scientific, California, USA), 5.25 µl of ddH 2 O and 10 ng of the template DNA. The cycling conditions were: one cycle at 95 °C for 5 min; 35 cycles at 95 °C for 30 s, 58 °C for 30 s and 72 °C for 30 s; followed by one cycle at 72 °C for 20 min. The PCR was performed in a Labnet™ MultiGene™ OptiMax Thermal Cycler (Labnet International, Inc.). The PCR product was purified by adding 0.25 µl of 10 U Exonuclease 1, and 2 µl of 2 U FastAP Thermosensitive Alkaline Phosphatase (ThermoFisher Scientific, California, USA) to the PCR product. The purification reaction was run for one cycle at 37 °C for 15 min followed by one cycle at 85 °C for 15 min in a Labnet™ MultiGene™ OptiMax Thermal Cycler (Labnet International, Inc.).
Cycle sequencing reactions were completed using the BigDye Terminator v3.1 Cycle Sequencing Kit (ThermoFisher Scientific, California, USA) using the Sanger chain termination method. Sequencing was conducted in a final reaction volume of 10 µl, containing 0.7 µl BigDye, 2.25 µl BigDye Terminator v3.1 5x Sequencing buffer, 0.75 µl ddH 2 O, 3.2 µM primer, and 5 µl of the PCR product. The cycling conditions were: one cycle at 95 °C for 2 min; 40 cycles at 85 °C for 10 s, 55 °C for 10 s and 60 °C for 2 min and 30 s. Reactions were performed in a Labnet™ MultiGene™ OptiMax Thermal Cycler (Labnet International, Inc.). Sequencing reactions were completed in both the forward and reverse direction. The cycle sequencing product was purified using the BigDye XTerminator Purification Kit (ThermoFisher Scientific, California, USA) and sequences were visualised using the 3500 Genetic Analyzer (ThermoFisher Scientific, California, USA).

Phylogenetic and genetic diversity analyses
Forward and reverse sequences were aligned in BioEdit (Hall 1999) to create a consensus sequence, and all sequences were manually trimmed and checked for ambiguous peaks. Absence of nuclear mtDNA sequences (numts) was confirmed following comparisons with published mitochondrial genomes and via visual inspection of the chromatograms. We included one ingroup reference sequence of H. senegalensis from Gabon (complete mtDNA genome; MN356338) and outgroup reference sequences from the genus. The three mitochondrial genes were initially analysed separately and included outgroup reference sequences of White-throated Kingfisher H. smyrnensis (KY940559; CO1, Cytb and 16S), Black-capped Kingfisher H. pileata (KJ476742; CO1 and Cytb) and Ruddy Kingfisher H. coromanda (MK327578 and KT356219; CO1, Cytb and 16S). For the concatenated dataset of the three mitochondrial genes (CO1, Ctyb and 16S), the outgroup reference sequences included a consensus of two species: the Black-capped Kingfisher H. pileata (KJ476742) and the White-breasted Kingfisher H. smyrnensis (KY940559). Inter-and intraspecific p-distances between subspecies and lineages for the concatenated mtDNA dataset were calculated using maximum likelihood (ML) genetic distance in MEGA7. Haplotype diversity (h), nucleotide diversity (π) and levels of gene-flow were all calculated using DnaSP software (Rozas et al. 2003), while Tajima's D statistic was estimated in Arlequin 3.5.2.2 (Schneider et al. 2000). The mismatch distribution (Rogers and Harpending 1992) of the H. senegalensis population was calculated using Arlequin 3.5.2.2 for both gene regions to test for demographic changes. The Harpending's raggedness index (HRI) (Harpending et al. 1993) and sum of squared deviations (SSD) were calculated to determine the fit of the expected frequencies under the demographic expansion model; a small value of HRI suggests an expanded population, while a large value indicates a stationary population (Harpending 1994). Analysis of molecular variance (AMOVA) was implemented in Arlequin 3.5.2.2 to estimate population differentiation by testing hypotheses about genetic variation between the two subspecies, and about geographic differentiation among the three localities.
The two nuclear genes were also analysed separately.  MG001221). The Grey-headed Kingfisher was the only available outgroup DNA sequence for the concatenated nuclear DNA genes (RAG1 and FIB5) partition (MG001221 and MG008244). All reference sequences were obtained from the National Center for Biotechnology Information (NCBI) sequence database GenBank. We determined the best-fitting substitution model using MEGA7 (Kumar et al. 2016). Haplotype reconstruction for the nDNA sequence was done after phasing the data in DnaSP, using the algorithms provided in Stephens et al. (2001) and Wang and Xu (2003) for heterozygous sites (polymorphic nucleotide positions).
We computed phylogenetic trees in MrBayes (Ronquist and Huelsenbeck 2003) and MEGA7, by using Bayesian inference (BI) of phylogenetic networks and ML methods for the separate as well as concatenated mtDNA and nDNA datasets. The BI analysis was conducted using the general time-reversible (GTR) model with the following parameters: the ML model employed six substitution types, with gamma-distributed rate variation (0.07) across sites, modelled using gamma partitioning of the mtDNA data only (mtDNA analysis and concatenated dataset of all genes). The Markov chain Monte Carlo (MCMC) approach was used to search for trees using four chains for 5 000 000 generations, with trees sampled every 50 000 generations (the first 20% were discarded as burn-in). The resultant trees were formatted in MEGAX (Kumar et al. 2018) with posterior probabilities. Phylogenetic relationships were also reconstructed using the ML method based on the Hasegawa-Kishino-Yano (HKY+G) model for CO1, Cytb and the concatenated mtDNA dataset; the Kimura 2-parameter (K2) model for 16S and RAG1; and the Jukes-Cantor one-parameter model for FIB5. Branch support values were estimated using nonparametric bootstrap with 1 000 replicates. The consistency of the two methods (BI and ML) was evaluated by comparing support values at identical nodes and by assessing the similarity of the tree topologies. Phylogenetic analysis of the concatenated mtDNA gene regions identified divergence between subspecies as well as demonstrated well-supported geographic lineages. The degree of divergence between subspecies and lineages was thus calculated using ML genetic distance in MEGA7.

Haplotype network inference
Haplotype networks were constructed separately for the concatenated mtDNA gene regions (CO1, Cytb, 16S) as well as for concatenated nDNA gene regions (RAG1, FIB5) using PopART v1.7 (Leigh and Bryant 2015), which employs an agglomerative approach where clusters are progressively combined with one or more connecting edges. This was done using the TCS haplotype network, a parsimony-based method of analysis defined by Templeton et al. (1992) and Clement et al. (2002).

Divergence-time estimation
Phylogenetic reconstruction and the calibration step in molecular dating analysis were done using Bayesian evolutionary analysis in BEAST v1.8.4 (Drummond et al. 2012) implemented with BEAGLE (Ayres et al. 2011) using the mtDNA datasets. A configuration file generated in BEAUti v1.8.4 was run using a random starting tree under the following parameters: GTR model, 10 million generations of the Markov chain, 10 000 tree-sampling frequency, with the first 10% discarded as burn-in, a strict clock model and speciation (Yule) process. The divergence times were calibrated using published date estimates for the three major clades of kingfishers of the family Alcedinidae (i.e. Alcedininae, Cerylinae and Halcyoninae). An average node age of 16.3 Ma (95% CI 13.2-19.6) was used for the root of the tree, representing the Miocene origin of the Halcyoninae (Anderson et al. 2018). Three species of its sister clade Cerylinae (i.e. Ceryle rudis, Megaceryle lugubris and Chloroceryle aenea) were included as the outgroup for this analysis. For the clade Halcyoninae, published mtDNA sequences of Halcyon (coromanda, pileata and smyrnensis) and haplotypes of H. s. senegalensis and H. s. cyanoleuca generated by this study were included in the analysis with three species from closely related genera (Todiramphus, Dacelo and Actenoides). (Information on the additional GenBank sequences used in this phylogenetic tree is given in Supplementary Table S1.) The information contained within the sampled trees of the BEAST output was summarised using TreeAnnotator v1.8.4, with the annotated tree visualized using FigTree v1.4.4 (http://tree.bio.ed.ac.uk/ software/figtree).

Molecular sexing
Molecular sexing was successful for 59 of the 60 birds tested. Thus, the total collected sample consisted of 30 males and 29 females: 17 females and 20 males from South Africa; 6 females and 4 males from Ghana; and 6 females and 6 males from Uganda.

Morphometric analyses
Morphometric measurements were collected for 52 adult individuals: 15 females and 20 males birds from South Africa, 6 females and 4 males from Ghana, and 4 females and 3 males from Uganda (Table 1a). Measurements taken from 5 juvenile birds were excluded. ANOVA results for effect of measurement errors were not significant (p > 0.05) for three variables (head length, tail length, and mass), showing no differences among measurement means taken by different samplers [F(2,45), critical value = 3.20, α = 0.05). However, there was a significant difference in measurements taken for tarsus length and wing length among samplers, but only owing to sampling in South Africa (Table 1a). Therefore, this difference was not attributed to measurement error but to actual differences in bird sizes among the regions. The PCA of the five variables recovered six factors that had eigenvalues of >1, explaining 99% of the variation for log-transformed and untransformed datasets (Table 1b). Principal component 1 (PC1) described 66.27% and 51.79% of the variation, while PC2 accounted for 19.86% and 26.04% of the variations for untransformed and log-transformed datasets, respectively. Morphometric variables with significant factor loadings in PC1, PC2 and PC3 were the mass, wing length and tail length, whereas mass, BMI and tarsus length remained significant after log-transformations. These variables all had factor loadings of more than 0.50 above the 0.30 that is considered significant. MANOVA indicated that there was no significant difference among sexes among localities for the two subspecies (Wilks lambda: 0.78, p = 0.07). A plot of PC1 and PC2 (Figure 2) showed that the male and female samples could not be segregated into groups, with a high degree of overlap. However, there was clustering among populations, with the Ghana individuals being more morphologically discrete from the other two populations (having little to no overlap), and with a higher degree of overlap between Uganda and South Africa.

Mitochondrial DNA phylogeographic and phylogenetic analysis
In this study, we generated sequences for 60 samples, of which 49 samples included all three gene regions (Supplementary Table S1): GenBank accession numbers OL602280-OL602343 for Cytb; GenBank accession numbers OL518993-OL519047 for CO1; and GenBank accession numbers OL519049-OL519104 for 16S. Lack of amplification of certain gene regions may be due to degradation of DNA at primer binding sites.

Analysis of subspecies
A total of 19 haplotypes were identified based on the concatenated mtDNA dataset, of which 9 haplotypes were unique to H. s. cyanoleuca (South Africa) and 10 were unique to H. s. senegalensis (Figure 3a). Singlelocus tree topologies showed two resolved monophyletic, divergent clades for all three mitochondrial gene regions (Supplementary Figure S1a Figure S2)    bootstrap support values (100% for both the ML and BI analysis). The overall dataset had very high haplotype diversity (0.73) and low nucleotide diversity (0.0007) ( Table 2). Tajima's D estimates were not significant (p > 0.05); however, the estimate was negative (−1.16) for H. s. cyanoleuca, suggesting that the South African population may have recently undergone a population expansion (Table 2). This result was also supported by the mismatch distribution (Supplementary Figure S4) for the species, which was bimodal and significant (p < 0.05) for the mtDNA data, suggesting either a recent population expansion or balancing selection favouring genetic variation among populations (Harpending et al. 1993). The nDNA had a more ragged mismatch (to the right, and broader) that was not significant (p > 0.05), displaying a high level of nucleotide variation, with individual sequences differing at many sites-a pattern expected following a bottleneck in a previously large population-with possible incomplete lineage sorting (Harpending et al. 1993). The low and non-significant estimates of HRI and SSD suggest the presence of non-equilibrium and that the nDNA data had a relatively good fit to the expansion model. Pairwise sequence divergence (F ST ) for concatenated mtDNA sequences between the two subspecies was very high and significant (>0.94, p < 0.05) (results not shown). However, the hierarchical AMOVA analyses of populations (locality groups) between the two subspecies indicated no significant variation (p > 0.05), although 83.88% of the total genetic variation was contributed by 'among groups' variation, with genetic variation from 'among populations within groups' accounting for 9.69% of the variation (Table 3). A significant 6.43% of the variation was detected within populations. Our analysis showed that the sequence divergence between H. smyrnensis and H. senegalensis, and between H. pileata and H. senegalensis, was 9.3% and 8.6%, respectively (    (5), Gabon (1), and Ghana (4) (Figure 3a; Table 2). Thus, no shared haplotypes were detected among the three H. s. senegalensis localities, suggesting geographic sub-structuring. The phylogenetic trees (CO1, Cytb and concatenated mtDNA) also showed sub-structuring in the H. s. senegalensis clade, with two lineages comprising samples from Uganda, and the second lineage including samples from Ghana and Gabon, with high bootstrap support values for the concatenated mtDNA (97% ML and 100% for BI analysis), suggesting isolation over the geographic range between the western and eastern regions of Africa (Figure 3c; Supplementary Figure S2). However, absence of sub-structuring was detected in the 16S phylogenetic tree (Supplementary Figure S1c). Pairwise sequence divergence (F ST ) for concatenated mtDNA sequences was high and significant (0.802, p < 0.05) between the two H. s. senegalensis geographic lineages ( Table 2). Distribution of the two lineages clearly corresponds to their geographic origins: Uganda and Ghana. AMOVA comparisons were done for the three groups, denoted 1 = South Africa, 2 = Uganda, and 3 = Ghana and Gabon, and further highlighted that there was higher genetic variation within populations than among assigned groups (Table 3). The estimate of gene flow (Nm) was 0.22, suggesting that the genetic differentiation among these populations is significant and that gene flow is restricted. Govindajuru (1989) indicated that levels of gene flow with Nm of <0.25 represent low gene flow, while Nm of >1 can be categorised as high gene flow, and 0.25-0.99 represents intermediate gene flow. Sequence divergence between geographic lineages within H. s. senegalensis (Ghana, Gabon and Uganda) varied between 0.2% to 0.4% (Table 4).

Nuclear DNA phylogenetic analysis
Here, we generated 704 bp of sequence for RAG1 (GenBank accession nos. OL602344-OL602402), and 801 bp of sequence for FIB5 (GenBank accession nos. Table S2b), represented by 92 phased sequences. Bayesian inference (BI) of the concatenated nuclear genes (1 505 bp) showed short internal branches, possibly caused by a rapid radiation of lineages. The BI tree (Supplementary Figure S2) and individual ML trees (Supplementary Figure  S1d,e) also showed a different topology to the mtDNA tree, and some lineages were not monophyletic. The overall data set has high haplotype diversity (0.99) and low nucleotide diversity (0.0113) ( Table 2). A total of 77 unique haplotypes for this dataset were identified, and none of the haplotypes were shared among the three population localities, similar to the mtDNA (Supplementary Figure S3) Table 4: Estimates of percentage sequence divergence between species, subspecies and geographic lineages of Halcyon kingfishers. SE estimates are shown above the diagonal 77 haplotypes that were not shared between subspecies or localities. Thus, biogeographic discordance and phylogenetic incongruence between mitochondrial and nuclear markers was observed in this study. Several hypotheses have been put forward to explain mitochondrial-nuclear discordance, including incomplete lineage sorting, hybridisation and ancestral population structure (Toews and Brelsford 2012;Linck et al. 2019). Funk and Omland (2003) reported that a high proportion of bird species (16.7%) were paraphyletic, with the most-common reason identified as incomplete lineage-sorting caused by recent speciation (McKay and Zink 2010  Periodic climate oscillations across the African continent have been reported, with three peaks of aridification being described at approximately 2.8, 1.7 and 1.0 MYA (deMenocal 1995). During these periods, repeated oscillations in average temperature and rainfall shifted between humidwarm (pluvial) and arid-cool (interpluvial) phases. During pluvials, forests expanded, and during interpluvials, arid areas expanded. In this study, we suggest that the common ancestor of these two subspecies was more broadly distributed in Africa ~1-2 MYA. The interpluvial period that occurred 0.8-1.2 MYA resulted in forest contraction and the expansion of savannah and grassland habitat (deMenocal and Bloemendal 1995), which is unfavourable to the Woodland Kingfisher. The pronounced dry period most likely closed the corridor between southern and eastern populations of Woodland Kingfishers, resulting in morphological/genetic divergence and differentiation. The reported arid period (0. . Here, we further suggest that contemporary restrictions to gene flow between Woodland Kingfisher subspecies may be due to allochrony, or divergence in breeding time (Servedio et al. 2011). This is associated with temporal variations, the availability and abundance of insects, geographic distances and/or local adaptation to the habitat. Woodland Kingfisher populations in South Africa are reported to breed in November, West African populations breed in June, and East African populations breed in January. Thus, seasonal separation of breeding times may be an important driver contributing to continued isolation of the subspecies (Taylor and Friesen 2017).

OL602403-OL602453) in 46 samples (Supplementary
Genetic structure was also recovered within H. s. senegalensis, suggesting that diversity within the Woodland Kingfisher is underestimated. Phylogenetic and network analyses of the mtDNA sequence data revealed significant differences at different geographic scales (Ghana, Gabon and Uganda), with three distinct lineages detected within this subspecies, suggesting a significant historic isolation among these populations. The separation between the lineages of H. s. senegalensis from western/ central Africa from eastern Africa occurred ~0.22-0.57 MYA, whereas the separation between western and central African populations occurred more recently at ~0.12-0.36 MYA. Divergence of populations from western, central and eastern Africa has been reported for several terrestrial vertebrates (lizards, mammals and birds) within a similar period. For example, chimpanzees constituted a single population until ~0.1 MYA, and were subsequently divided into three populations (southern Cameroon, central Africa, and eastern Africa) (Gonder et al. 2011). Two species of woodpeckers (Campethera caroli and C. nivosa) from the upper and lower Guinean forest blocks were reported to each consist of populations that diverged between 0.5 and 0.8 MYA (Fuchs and Bowie 2015). Furthermore, Perktaş et al. (2020) described significant genetic and morphological variation within turacos (family Musophagidae) that were generally concordant with the organisation of the montane avifaunal regions of Africa, which have in the past been driven by Pliocene Forest dynamics. Several potential biogeographic barriers have been identified, including the Dahomey Gap, which has been described as a 200-km wide forestsavannah mosaic corridor separating the rainforests of West Africa and central Africa (Moreau 1966;Hall and Moreau 1970;Crowe and Crowe 1982). The habitat and size of the Dahomey Gap has varied over time in response to large-scale climate shifts (Mayr and O'Hara 1986). Additional potential barriers include the Niger Delta and the Cameroon volcanic line, which are reported to be barriers to gene flow between populations from West Africa and Central Africa (Hassanin et al. 2015). Lastly, the Congo rainforest may potentially present a source of isolation between the central (Gabon) and eastern African group.
The morphological analysis supported the structuring results detected with mtDNA, and our results showed morphological differentiation between the eastern, western and southern populations. Limited studies to date have reported on morphological variation within Woodland Kingfisher populations and have only described phenotypic differences between subspecies, such as the presence of an eye stripe in H. s. cyanoleuca (Fry et al. 1988). In our study, birds from South Africa had longer wing and tail lengths, had lower mass, and constituted a highly differentiated genetic cluster. The lower mass observed in H. s. cyanoleuca may be attributable to physiological processes that occur as a result of long-distance migration (Ramenosky 1990) and that can vary during an annual cycle. Birds from Uganda were intermediate in those measurements, whereas Woodland Kingfishers from Ghana had the shortest wing and tail lengths. The H. s. senegalensis in southern parts of Ghana are reported to be resident populations; hence, shorter wings in populations that do not migrate is expected (Perez-Tris and Tellerıa 2001). However, in this study we additionally detected a continental gradient in wing length, from the western to the eastern populations.
The molecular analysis confirmed the subspecies designations of H. s. cyanoleuca and H. s. senegalensis. In addition, we further identified distinct lineages within H. s. senegalensis. Here, we provide evidence that climate change leading to expansion and contraction of geographic range has played an important role in shaping populations of Woodland Kingfishers. We further suggest there is limited contemporary gene flow in these populations caused by distance, differences in breeding behaviour, and/or local adaptation to their habitat. Our results concur with findings for other vertebrate species that have identifiably unique populations in central, eastern, western, and southern Africa. However, additional data need to be collected throughout the distribution range of Woodland Kingfishers to develop a comprehensive picture of intra-and interspecific variation within H. s. senegalensis and how this is distributed geographically with the other co-occurring subspecies. The relationships of H. senegalensis, including the H. s. fuscopileus subspecies that we were not able to sample in our study, clearly awaits molecular analysis. Although Woodland Kingfishers are currently abundant and have been assessed as Least Concern, the identification of unique populations in a continually transforming habitat may require future conservation management.