Boom et al. (2021) provided the first phylogeny of the genus Brachystegia by sequencing full plastomes but plastid haplotypes appeared shared among species at a local scale, suggesting widespread plastid capture (and thus, hybridization to some extent). Consequently, the plastid phylogeny reflects geographical rather than taxonomical affinities between samples, a pattern well documented for Quercus, Fraxinus, and Macaranga trees (e.g. Bänfer et al. 2006; Heuertz et al. 2006; Simeone et al. 2018). The plastid phylogeny delineates two main parapatric clades separating most of the Guineo-Congolian specimens from the Zambezian ones (Boom et al. 2021). The Zambezian clade is structured in three parapatric clades, ranging from Tanzania to Angola (with eastern, central, and western subclades). The first cladogenesis event occurred during the late Miocene-Pliocene, while the different Zambezian plastid clades originated in the Pliocene-Pleistocene. Interestingly, a longitudinal gradient of time to the most recent common ancestor (TRMCA) is observed for the three Zambezian subclades and suggests a westward expansion of the Miombo woodlands, from an original range situated in East Africa. Although plastid clades have provided many insights into the biogeographic history of Brachystegia, the evolutionary relationships between the species remain to be characterised. To this end, phylogenetic information from the nuclear genome is required.

The nuclear-encoded 35S ribosomal DNA (rDNA) cistron, comprising the 18S, 5.8S, and 25S rRNA genes, is the most easily accessible nuclear DNA region using genome skimming, i.e. shallow genomic sequencing (Straub et al. 2012). However, interpreting variation in 35S rDNA sequences and their internal transcribed spacers (ITS1 and ITS2) can be challenging due to their mode of evolution (Feliner and Rosselló 2007). In plants, the nuclear-encoded 35S rDNA is located in one or several loci where hundreds of copies follow each other in arrays that tend to be homogeneous at the intra-individual scale due to concerted evolution (Eickbush and Eickbush 2007). Complete concerted evolution is however not universal (Bailey 2003). Hence, a certain degree of polymorphism can occur inside or between arrays, due to heterozygosity, homeology, or paralogy, including pseudogenes (array silencing, e.g. Volkov et al. 2017), causing Intra-Individual Sites Polymorphisms (2ISPs; Potts et al. 2014). When 2ISPs are not due to heterozygosity, they can sometime be fixed within a species, and therefore be phylogenetically informative. A potential trade-off between accuracy and simplicity in data analysis is to encode 2ISPs using IUPAC codes and to evaluate their putative impact on phylogenetic reconstructions. Potts et al. (2014) showed that treating 2ISPs as informative states rather than as ambiguous or missing characters could increase the resolution of phylogenetic reconstructions for highly polymorphic datasets (Potts et al. 2014, but see Fonseca and Lohmann 2020).

In addition to 2ISPs, other processes could blur the phylogenetic relationships between species at the nuclear genome, namely hybridization and incomplete lineage sorting (ILS) (Maddison 1997) or lack of divergence (e.g. Turner et al. 2016). Hybridization is suspected to occur among Brachystegia savanna species based on morphological arguments (White 1962; Brenan 1967) but also given the evidence of recurrent plastid capture among species (Boom et al. 2021). ILS could also be likely, as dominant savanna tree species may have diversified recently and/or have large population sizes and long generation times, resulting in retention of ancestral polymorphism (Pennington and Lavin 2016). Finally, when DNA sequences are too short and/or evolve too slowly to accumulate enough mutations along the branches of a species tree, they offer limited resolution by lack of divergence.

This work aims to investigate the diversification of Brachystegia, bridging the gap between genetic and morphological/taxonomic considerations through a multimarker phylogenetic approach. Namely, we use genome skimming on specimens from museum collections with degraded DNA (Zeng et al. 2018; Alsos et al. 2020) to assemble their nuclear ribosomal DNA, supplementing insights previously obtained from their plastomes (Boom et al. 2021). We will address the three following questions. Firstly, is the rDNA region phylogenetically informative for reconstructing the Brachystegia phylogeny? Secondly, does taking in account 2ISPs variability increase the resolution of our phylogenetic reconstructions? And thirdly, does rDNA provide complementary information to plastid DNA to infer the evolutionary history of the genus?

Brachystegia plant material (leaves) was collected on vouchers from the four following herbaria: BR, BRLU, FHO, and LISC (acronyms according to Index Herbariorum, Thiers continuously updated). In addition, we collected material (n = 10) from fieldwork and dried leaves using silica gel for DNA extraction purposes. In total, 47 individuals of 27 of the 29 described species of Brachystegia (following Lebrun and Stork 2008) were sequenced (Table 1). A sample of Julbernardia paniculata was added as outgroup to root the different reconstructed phylogenies. On both silica-dried and herbarium leaves, DNA extractions were performed using the DNeasy Plant Mini Kit (Qiagen, Netherlands) and the protocol detailed in Cappellini et al. (2010), except that the digestion step was performed overnight at 37°C rather than 55°C and we did not perform an initial wash step of the plant material with a bleach solution. DNA concentration and DNA size distribution were assessed with a Qubit® 2.0 Fluorometer (Thermo Fisher Scientific, USA) and with an electrophoresis on a 1% agarose gel. The different genomic libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs, USA) and were pooled equimolarly. After pooling, we sequenced the libraries on an Illumina NextSeq 500 instrument at the GIGA platform (Liège, Belgium), using the V2 mid-output reagent kit and targeting one million reads per library (2 × 150 paired-end reads).

First, the quality of the genomic libraries was checked using FastQC v.0.11 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and then they were trimmed to remove low-quality regions and Illumina adapters using Trim Galore! v.0.4.5 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). De novo assemblies of the rDNA sequences were done for each sample using GetOrganelle v.1.6.2 (default parameters for ribosomal assembly; Jin et al. 2020). Resulting graphs were inspected using Bandage v.0.8.1 (Wick et al. 2015). The rDNA sequences were annotated using Infernal cmscan (Madeira et al. 2019) to delineate 18S rDNA, internal transcribed spacer 1 (ITS1), 5.8S rDNA, internal transcribed spacer 2 (ITS2), and 25S rDNA regions of the 35S cistron. Reads were finally mapped on their respective assembly using the Burrows-Wheeler Aligner BWA mem v.0.7.12 (Li and Durbin 2009), and depth coverage was computed using Samtools v.1.9 (Li et al. 2009). We polished the assembled rDNA sequences and identified 2ISPs positions by conducting a consensus calling for each individual with the following parameters: min-var-freq 0.01, --min-avg-qual 30, --min-freq-for-hom 0.9, and --min-coverage 3, using Varscan v.2.3.7 (Koboldt et al. 2012). Polymorphic sites were coded using IUPAC recommendations. Additionally, indels were coded as N. The polished sequences were deposited in GenBank (Table 1). The 18S–25S rDNA were aligned using MAFFT v.7 (Katoh et al. 2019) and the alignment was visually checked using MEGA7 (Kumar et al. 2016). The alignment is available in Supplementary file 1.

In order to evaluate the impact of 2ISPs, we performed Maximum Likelihood (ML) based phylogenetic inferences with RAxML-NG software using three different implementations (Kozlov et al. 2019). (i) All 2ISPs were recoded as missing data (ML-N for missing data ‘N’). (ii) Binary 2ISPs were treated as ambiguous with both states of the coded ambiguity considered equiprobable, while 2ISPs involving three nucleotides were coded as missing data, using a GTR+I+G substitution model (ML-A for ‘Ambiguous’; Potts et al. 2014). (iii) Finally, we considered each 2ISP genotype as a transitional state, with the transitional probability included as additional parameter into the substitution model. This approach (ML-I for ‘Informative’) followed a comparable logic as in Potts et al. (2014), except that we preferred the new GENOTYPE implementation in RAxML-NG, using the GTGTR4+I+G substitution model (see Kozlov et al. 2022: supplementary note 2). For each method, branch supports were computed based on 1000 bootstrap replicates and are provided as nonparametric bootstrap support (BS).

Mutation patterns and rates differ between the rRNA genes and the internal transcribed spacers (Baldwin et al. 1995). To assess the robustness of our phylogenetic reconstructions, we also used RAxML-NG with the different implementations described above but using two partitions (ITS1+ITS2 vs all rRNA genes; GTR+I+G model for each partition in ML-N and ML-A; GTGTR4+I+G model for ML-I). We also evaluated how the different main subunits of the 35S cistron are phylogenetically informative and/or might provide conflicting signals by conducting separate RAxML-NG analyses on each region (ITS1, ITS2, 18S, and 25S) with the GTGTR4+I+G substitution model. Finally, at shallow phylogenetic levels, we can expect non-tree-like evolutionary patterns (i.e. not explained by strictly bifurcating trees) due to reticulation, low-level genetic divergence and ILS. We therefore also applied network approaches and computed a Neighbor-Net splits graph using p-distance (i.e. 2ISPS aware approach) computed with the R package phangorn v.2.5.5 (Schliep 2011; Schliep et al. 2017) and SplitsTree5 (Huson and Bryant 2006). We also computed a bootstrap consensus network (e.g. Schliep et al. 2017) using the ML-I bootstraps trees (ML-I 1 partition, 1000 BS trees; edge weights: tree size weighted mean; Threshold: 15).

Finally, to evaluate the ability of rDNA to retrieve clades with taxonomic significance, we assessed whether sequences divergence was correlated to taxonomy and/or geographical distance using Mantel tests computed with the R package vegan v.2.5-6 (Oksanen et al. 2019). For the different clades identified in the reconstructed phylogenies, we searched for correlation between matrices of genetic distances, computed with the R package ape v.5.4-1 (raw distance; Paradis and Schliep 2019), and taxonomic distances (0 and 1 indicating whether two individuals belong to the same species or not). Additionally, as geographic and genetic distance correlation can provide information regarding the geography of diversification (Abellàn and Ribera 2017) and/or events of hybridization at a local scale (e.g. Boom et al. 2021), Mantel tests were also performed between geographic (i.e. shortest distances computed by the R package geosphere v.1.5-10; Hijmans 2019) and genetic distances (R Core Team 2019).

After the trimming step, we obtained a mean of 1,227,134 reads per library (SD = 642,412). Raw sizes of the different rDNA assemblies ranged between 5,948 bp and 9,048 bp (mean = 7,526, SD = 757), and the length of the extracted 18S–ITS1–5.8S–ITS2–25S sequences ranged between 5,823 bp and 5,843 bp (mean = 5,828 bp). Mean depth coverage per sample ranged between 49X and 1,096X (mean = 260X, median = 192X). This wide range in mean depth coverage was also reflected in terms of relative rDNA reads content, i.e. between 0.19% and 3.35% of the total number of reads for each genomic library mapped on the corresponding assembly (mean = 0.81, SD = 0.66). The final alignment included a total of 5,973 sites and contained 166 variable and 78 parsimony-informative sites, when 2ISPs were not taken into account. When the outgroup J. paniculata was removed, these values drop to 108 and 75 variable and parsimony-informative sites, respectively. Taking into account the exclusive 2ISP sites, i.e. sites characterised by a single-defined nucleotide with ambiguity codes, the number of variable sites rose to 394 and 338, respectively with and without J. paniculata. Details regarding the exact base composition of each sequence are available in Supplementary file 2.

ML-N and ML-A phylogenetic reconstructions supported similar groupings, with minor differences (more details in Supplementary file 3: Fig. S1). For clarity, we only present the ML-A and ML-I cladograms, i.e. without branch lengths, in the main article (Fig. 1). The ML-N, ML-A, and ML-I phylograms are available in Supplementary file 3 (Figs S2, S3, and S4). Both ML-A and ML-I inferences provided topologies with low bootstrap (BS) supports for most branches (Fig. 1). Additionally, among the 16 species represented by at least two individuals, only for four (ML-A) or six (ML-I) of those species, the different individuals were placed together. At the intra-generic level, we delineated four robustly supported clades that form coherent taxonomic, ecological, and/or spatial entities. One to two clades of rain forest species (B. eurycoma-B. nigerica and B. kennedyi-B. leonensis clades) and two clades of Miombo species are well supported (BS = 75–100 for ML-A, 84–99 for ML-I; Fig. 1). The positions of three species, B. cynometroides, B. laurentii, and B. mildbraedii, remain unresolved, even if the ML-I analysis suggested a B. cynometroides-B. laurentii grouping (BS = 64). The two identified Miombo clades are morphologically distinct according to floral, leaf, and axillary dormant bud features, corresponding roughly to the morphogroups A and B defined in the Flora of Tropical East Africa (Brenan 1967). The Miombo Group A clade includes B. bakeriana, B. bussei, B. floribunda, B. manga, B. microphylla, B. puberula, B. spiciformis, B. tamarindoides, B. torrei, and B. utilis, whereas the Miombo Group B clade includes B. allenii, B. angustistipulata, B. boehmii, B. gossweileri, B. longifolia, B. russelliae, B. stipulata, B. taxifolia, and B. wangermeeana. Except for B. russelliae in both the ML-A and ML-I method and for B. boehmii in ML-I, none of the Miombo species was sorted taxonomically. In addition, three species with similar leaf morphology (i.e. B. bakeriana, B. floribunda, and B. spiciformis) are placed together but with limited support in the three different implementations (BS = 40 for ML-N, 44 for ML-A, 39 for ML-I).

Partitioning the data (ITS1+ITS2 and rDNA genes) did not produce major differences that are supported with high bootstrap values (Supplementary file 3: Figs S5, S6, and S7b). We again observed three to four robustly supported clades, low BS values for most of the branches, and individuals from the same species were only placed together in four or six cases depending on the RAxML-NG implementation used (Supplementary file 3: Figs S5, S6, and S7).

The Neighbor-Net splits (Fig. 2) and the bootstrap consensus network graph (Supplementary file 3: Fig. S8) also identified four main clusters corresponding to the four main clades (i.e. Miombo Group A, Miombo Group B, B. eurycoma-B. nigerica, and B. kennedyi-B. leonensis). ITS2 was the most informative subunit as it identified the four main clades, contrary to the other subunits, although the bootstrap support values were moderate to low (BS = 19–65; Fig. 2). The least informative subunit was the 18S rDNA gene, where B. eurycoma-B. nigerica is the only main clade found, with moderate support (BS = 60; Fig. 2). The individual trees are available in Supplementary file 3 (Figs S9–S12).

The correlation between the genetic and taxonomic distances was significant for the Miombo species as a whole (r = 0.19, p value = 0.001, n = 35) but not within Miombo Group A (r = 0.11, p value = 0.138, n = 17) or Miombo Group B (r = -0.069, p value = 0.219, n = 18). Hence, rDNA discriminated between different Miombo taxonomic groups but did not provide fine-scale taxonomic information. The correlation between geographical and genetic distances was non-significant for Miombo species as a whole (r = 0 .12, p value = 0.028, n = 35) and for Miombo Group B (r = -0.02, p value = 0.551, n = 18) but was significant within Miombo Group A (r = 0.34, p value = 0.001, n = 17). The mean genetic distance between individuals was higher in Miombo Group A than in Miombo Group B (i.e. differences in the nucleotide composition of 0.03% and 0.0055%, respectively) for roughly the same number of species (n = 10) and specimens (n = 17–18). Overall, the number of parsimony-informative sites within Miombo Group A and B was low (n = 17 in both clades).

Explicitly taking into account the retention of polymorphism allows for the detection of clades and species groups that align with morphology. For instance, regardless of the partitioning scheme (one vs two partitions), ML-N placed the specimen B. nigerica 30 (Lapido 19061, Nigeria) as poorly supported sister to the B. kennedyi-B. leonensis clade (BS = 22–28), while the specimen B. nigerica 29 (Chesters A124/30, Nigeria) grouped with B. eurycoma (BS = 61–63). In ML-A, the specimen B. nigerica 30 moved to the B. eurycoma-B. nigerica subtree but with decreased support (BS = 9). With two partitions, however, the specimen B. nigerica 30 remained sister, with low support, to the B. kennedyi-B. leonensis clade (BS = 25). In contrast, in ML-I, regardless of the partitioning scheme, both specimens of B. nigerica were resolved as sisters with moderate support (BS = 53–57) in a well-supported clade composed of B. eurycoma and B. nigerica (BS = 80–84).

As described in Boom et al. (2021), the plastid phylogeny delineates clades that correspond to geographic regions (Fig. 3). Five main plastid clades and two additional singleton lineages delineated seven regions. Two clades are part of a larger rain forest (RF) clade (Fig. 3; pentagon = Upper Guinea and Southwest Nigeria region; diamond = Southeast Nigeria-Cameroon region). The RF sister clade encompasses two basal lineages (Fig. 3; triangle = the Eastern Arc Mountains and surroundings; reverse triangle = Lower Guinea) in addition to three parapatric Miombo woodlands (MW) clades (Fig. 3; circle = Eastern; star = Central; square = Western). Brachystegia specimens in Miombo regions share plastid sequences from different clades, regardless of the rDNA clades to which they belong (most specimens from rDNA Miombo A and B clades possess plastid sequences from one of the three parapatric Miombo clades). Nuclear ribosomal DNA is therefore sorted according to the taxonomy (to some extent), while plastid DNA is geographically structured. This supports cytoplasmic genome exchanges between all the different Miombo species, including between Miombo Groups, as well as between RF species. For example, the specimen B. nigerica 30 has a plastid sequence that is very similar to those in the geographically nearby specimens B. kennedyi 17 and 18 (Upper Guinea and Southwest Nigeria clade), while B. nigerica 29 plastome is part of the Southeast Nigeria-Cameroon region plastid clade including the two geographically nearby B. eurycoma specimens (see Table 1). In contrast, the ML-I rDNA phylogeny supports the specimen B. nigerica 30 as being part of the B. nigerica-B. eurycoma clade, suggesting at least one event of B. kennedyi plastid capture by B. nigerica in Southwest Nigeria.

The sorting of individuals according to their taxonomic species in the rDNA trees is, for the Miombo woodland species, the exception irrespective of the treatment of 2ISPs. However, reciprocal monophyly is supported for several Guineo-Congolian species and for the morphologically distinct Miombo Group A vs Miombo Group B. Correlation between geographic and genetic distances for Miombo Group A could suggest some inter-species gene flow at local scale but could not rule out past allopatric diversification. On the other hand, such correlation was not found in the Miombo Group B. Altogether, rDNA sequences provide insights regarding the Brachystegia evolutionary history, even if most of the relationships between species and clades remained unresolved due to insufficient phylogenetic information. This lack of discriminating signal between species could be due to the fairly recent origin of the different species (Boom et al. 2021), which could have evolved in parallel with the origin and expansion of C4 fire prone savannah during the late Miocene-Pliocene-Pleistocene (Maurin et al. 2014; Polissar et al. 2019). The relationships between the four main clades of Brachystegia and several of the Guineo-Congolian species (e.g. B. cynometroides, B. mildbraedii, B. laurentii) remain unresolved here, possibly because these lineages diverged in a short time (i.e. rapid radiation).

Nuclear ribosomal DNA also retrieved clades that are congruent with the taxonomy. The formerly recognized Miombo A and B groups in F.T.E.A. (Brenan 1967) (hereinafter “morphological Group A and B”) seem to represent two monophyletic groups, consistently supported as clades in our different reconstructions. Being established for East African species, the infrageneric system of Brenan (1967) does not include per se all Zambezian woodland species. Namely, four western Miombo species (B. bakeriana, B. gossweileri, B. tamarindoides, B. russelliae) and three narrowly distributed species (B. michelmorei, B. oblonga, B. torrei) were not covered by Brenan (1967). However, these seven species were covered by the Flora Zambesiaca (Brummitt et al. 2007). They show clear morphological affinities with species included in Brenan (1967). Brachystegia tamaridoides and B. torrei share traits with B. microphylla, while B. bakeriana has many traits in common with B. spiciformis. On the other hand, B. gossweileri, B. russelliae, and B. michelmorei are, based on buds/flowers/stipules/auricles, part of the morphological Miombo Group B. Miombo Group species exhibit different vegetative and floral morphological trends (Brenan 1967). Morphological Group A species have globoid or ovoid buds, while species from morphological Group B have flattened buds enclosed in two large keeled scales. Morphological Group B species mostly have persistent stipules with basal reniform auricles, while in morphological Group A species stipules are generally caducous, mostly without auricles. Morphological differences are also reported for the bark, with relatively thick bark with vertical furrows in morphological Group B species. Other vegetative traits that have taxonomic value (e.g. the number and dimension of leaflets pairs) do not appear to show consistent differences between the two morphological groups. Regarding floral traits, tepals are absent, reduced, or shortly ciliate in morphological Group A species vs densely and long-ciliate in morphological Group B species. Most Brachystegia species have paniculate inflorescence, with the notable exception of three Miombo species: B. stipulata (morphological Group B) has raceme and/or paniculate inflorescences, while B. bakeriana and B. spiciformis species (morphological Group A) have racemose inflorescences (Brummitt et al. 2007). The latter two species are not always easy to discriminate in herbarium material (Brummitt et al. 2007), even if they differ in their habit, fruit, and bark traits (B. bakeriana is a shrub, while B. spiciformis is a tall tree) (White 1962; Brummitt et al. 2007). Sterile material of these species can be difficult to separate from B. floribunda (Brummitt et al. 2007). In the different implementations to reconstruct an rDNA tree (ML-A, ML-N, and ML-I), we systematically retrieved a moderate/low supported clade that includes individuals from the three aforementioned species. Such clade can reflect a close evolutionary relationship between these species. It is somehow unexpected, as B. floribunda is morphologically closer to other species from the Miombo Group A (i.e. several sepaloid/shortly ciliates tepals and panicles for B. floribunda and other Miombo group A species, while tepals are reduced or absent in other species). Alternatively, we cannot exclude that some of the specimens were not correctly identified and might contribute to this apparent and unexpected species cluster. Preliminary results using additional nuclear markers and a denser sampling suggests that the identification of B. floribunda specimens might be problematic (Boom 2021).

Apart from Zambezian species, the relationship between B. leonensis and B. kennedyi is congruent with the views of Hoyle (1955), who recognized close affinities of leaf anatomy, with however substantial differences in floral traits.

The rDNA phylogeny (Fig. 2) is in sharp contrast with the recently published plastid gene tree showing a strict geographical sorting (Boom et al. 2021). The large number of plastid introgression events observed in the genus resulting in large spatial clusters supports at least some interspecific gene flow within and even between the major groups. Hybrids are suspected to occur within each Miombo woodland group, while morphological intermediates between the two groups were rarely observed (Brenan 1967), in agreement with the rDNA data. Correlation between genetic and spatial distances was not found for Miombo Group B, and uncertainty on the exact reason for correlation in Miombo Group A prevents us to formally designate hybridization as the main driver of rDNA genetic diversity distribution among species. This may reflect allopatric speciation, hybridization, or even unbalanced sampling. However, low within-group resolution prevents us from pinpointing nuclear introgression between species.

Properly evaluating the effect of coding intra-individual polymorphisms is out of the scope of this paper. However, we note that considering the polymorphism as recommended in Potts et al. (2014) allowed to identify clades that are congruent with morphological observations and increase support of critical branches. Therefore, 2ISPs could contain valuable phylogenetic information, as in recognizing the B. eurycoma and B. nigerica grouping, which is significant given their frequency in rDNA sequences (we found about three times more variable sites among Brachystegia samples when including 2ISPs). The exact mechanisms producing 2ISPs remain to be characterized. Here, it is important to highlight that we mainly focused on tree topologies, and not on the other properties of phylogenetic trees, i.e. branch lengths. Works on Scrophularia and Bignoniaceae highlighted the strong impact of ambiguities coding schemes on branch lengths (Scheunert and Heubl 2017; Fonseca and Lohmann 2019) but they did not apply the here used RAxML-NG ML-I implementation. As seen in Supplementary file 3: Figs S2–S4, with our data, the treatment of 2ISPs did not substantially affect branch lengths despite a substantial impact on topology.

Overall, the rDNA gene trees in this study shed some light onto part of the evolutionary history of Brachystegia species but they cannot resolve relationships between closely related species because the assembled 18S–25S rDNA data is not sufficiently divergent and thus not sufficiently informative. The lack of resolution between Miombo species could be explained mainly by the non-mutually exclusive following reasons: a recent evolutionary history of diversification and/or the occurrence of gene flow between the extant species.

If due to recent diversification, the proper characterisation of a species tree may need the use of several unlinked loci. Targeted enrichment could constitute an interesting strategy, as such methods allow the investigation of young and species rich genera (e.g. Inga in Nicholls et al. 2015). Moreover, specific baits already exist for the Detarioideae subfamily and already proved their potential to unravel species trees in recent genera, allowing a proper characterization of both evolutionary history and taxonomy (Ojeda et al. 2019; de la Estrella et al. 2020). Alternatively, Genotyping by Sequencing (GBS) or restriction-site associated DNA sequencing (RAD-seq) have been proven useful in order to evaluate species trees for plants, even in presence of polyploidy, hybridization, and incomplete lineage sorting (e.g. Afzelia in Donkpegan et al. 2020; Cycnoches in Pérez-Escobar et al. 2020; Quercus in Hipp et al. 2020). Errors of identification in the diagnosed specimens (e.g. confusion between B. bakeriana, B. floribunda, and B. spiciformis) could be identified by such approaches, in addition by using a denser sampling for the problematic species (e.g. Inga in Dexter et al. 2010).

If the lack of resolution is due to reticulate history, the characterization of the nature, directionality, and extent of gene flow could be explored using targeted enrichment (e.g. Brownea in Schley et al. 2020). Moreover, the presence of interspecific gene flow could trigger a wider reflection on the species delineation within Brachystegia. The particularly weak nuclear genetic distances observed among the Miombo Group B individuals could potentially be explained by over-taxonomisation. Several species from the Miombo Group B form a morphological closely-related series with many intermediates (i.e. B. allenii, B. angustistipulata, B. boehmii, B. longifolia, and B. wangermeeana in Brenan 1967). Such a morphological continuum with different morphotypes might reflect shallow divergence, and the current taxa could be considered as having an infraspecific rank instead. Alternatively, the divergence between the current taxa could still be meaningful and all species would then be part of a wider interbreeding system (with limited gene flow between the different taxa), i.e. a syngameon (Grant 1981; e.g. oaks in Cannon and Petit 2020). The characterization of such a system could rely on combined genetic and morphological approaches (e.g. Tovar-Sánchez and Oyama 2004).