{'fraction_non_alphanumeric': 0.03450611408069122, 'fraction_numerical': 0.02460276781137878, 'mean_word_length': 5.117469204927212, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Between a Rock and a Hard Polytomy: Phylogenomics of the Rock-Dwelling Mbuna Cichlids of Lake Malaŵi

Abstract Whole genome sequences are beginning to revolutionize our understanding of phylogenetic relationships. Yet, even whole genome sequences can fail to resolve the evolutionary history of the most rapidly radiating lineages, where incomplete lineage sorting, standing genetic variation, introgression, and other factors obscure the phylogenetic history of the group. To overcome such challenges, one emerging strategy is to integrate results across different methods. Most such approaches have been implemented on reduced representation genomic data sets, but whole genomes should provide the maximum possible evidence approach. Here, we test the ability of single nucleotide polymorphisms extracted from whole genome resequencing data, implemented in an integrative genomic approach, to resolve key nodes in the phylogeny of the mbuna, rock-dwelling cichlid fishes of Lake Malaŵi, which epitomize the phylogenetic intractability that often accompanies explosive lineage diversification. This monophyletic radiation has diversified at an unparalleled rate into several hundred species in less than 2 million years. Using an array of phylogenomic methods, we consistently recovered four major clades of mbuna, but a large basal polytomy among them. Although introgression between clades apparently contributed to the challenge of phylogenetic reconstruction, reduction of the data set to nonintrogressed sites still did not help to resolve the basal polytomy. On the other hand, relationships among six congeneric species pairs were resolved without ambiguity, even in one case where existing data led us to predict that resolution would be difficult. We conclude that the bursts of diversification at the earliest stages of the mbuna radiation may be phylogenetically unresolvable, but other regions of the tree are phylogenetically clearly supported. Integration of multiple phylogenomic approaches will continue to increase confidence in relationships inferred from these and other whole-genome data sets. [Incomplete lineage sorting; introgression; linkage disequilibrium; multispecies coalescence; rapid radiation; soft polytomy.]

Complete genomes offer the tantalizing promise of resolving even the most recalcitrant phylogenetic relationships. However, rapid radiations with extensive incomplete lineage sorting (ILS) pose challenges due to their high gene-tree discordance (Edwards 2009). This results in unresolved or poorly resolved nodes in species trees. This is commonplace, and reconstructions of old radiations almost invariably contain one or more nodes that defy resolution, for example, in all major tetrapod groups (Irisarri and Meyer 2016;Suh 2016;Irisarri et al. 2017; Moreira and Schrago 2018;Braun et al. 2019;Hime et al. 2021;Singhal et al. 2021), and various plant groups (Pease et al. 2018;Smith et al. 2020;Cai et al. 2021;Gagnon et al. 2021;Morales-Briones et al. 2021). Although these problems often plague ancient radiations where they can be compounded by other effects such as saturation (e.g., Morales-Briones et al. 2021), polytomies can be just as pervasive in young, explosive radiations, and the problem is likely exacerbated when there is the opportunity for extensive gene flow (Koblmüller et al. 2010;Malinsky et al. 2018).
The cichlid fish radiations of the African Great Lakes epitomize recent and explosive adaptive radiations with gene-flow. These radiations have generated thousands of species in less than 10 million years (Meyer et al. 1990;Meyer 1993;Kornfield and Smith 2000;Seehausen 2006;Brawand et al. 2014;Henning and Meyer 2014;Malinsky et al. 2018) and have a high potential for hybridization (Kocher et al. 1995;Hulsey et al. 2010;Mims et al. 2010;Brawand et al. 2014;Genner et al. 2015;Svardal et al. 2021). As a result of their rapid diversification and frequent opportunity for introgressive hybridization, it is unclear whether the genomes of these fishes contain the information needed to parse their evolutionary history. It may even be impossible to recover a robust phylogenetic hypothesis for much of the diversity of these fishes. Alternatively, it might be that the relationships among the numerous genera and species could readily be placed within a strictly bifurcating and highly supported topology given sufficient genetic data and appropriate phylogenomic methods that take into account the various processes that obfuscate the reconstruction of their evolutionary history. Within this diversity of cichlid fishes, one of the greatest remaining phylogenetic challenges is found in the Lake Malaŵi rock-dwelling cichlids, the mbuna. These consist of 122 currently named species placed in 14 genera (Fricke et al. 2020) and are estimated to include hundreds of undescribed species that all arose in the last 2 million years (Meyer et al. 1990;Kocher et al. 1995;Genner et al. 2015;Schedel et al. 2019).
Previous attempts to reconstruct relationships among Malaŵi cichlids using high-throughput reduced representation sequencing techniques have had limited success. For instance, ultraconserved elements (UCEs) that can reliably isolate the same several hundred 741 742 SYSTEMATIC BIOLOGY VOL. 71 Notes: Clade is assigned according to our phylogenomic results, below. We here retain the species assignment of the taxa from the original studies, but note that several may have since undergone taxonomic changes, for example, Labeotropheus species.
loci produced incongruent phylogenies among different studies and phylogenetic reconstruction methods (McGee et al. 2016;Hulsey et al. 2017;Irisarri et al. 2017). With a chromosome-level assembly of the Malaŵi mbuna cichlid Maylandia zebra, coupled with the ability to resequence genomes with short reads, genome-wide information has now been devoted to the resolution of phylogenetic relationships among Malaŵi cichlids (Brawand et al. 2014;Malinsky et al. 2018;Conte et al. 2019;Urban et al. 2021). Yet, even with these nearly complete genomic data sets, incongruities still remain at several levels. The power to resolve relationships in an adaptive radiation that has undergone explosive diversification may depend on the stage of the radiation, or rough taxonomic level, at which relationships are examined (Albertson et al. 1999;Streelman and Danley 2003). Many of the mbuna genera and species are readily diagnosable phenotypically (Reinthal and Meyer 1997;Streelman and Danley 2003;Hulsey et al. 2010;Pauers 2010;Hulsey et al. 2013;Kratochwil et al. 2018). Therefore, wholegenome sequencing might be expected to be able to resolve species-and genus-level relationships among much of this radiation despite its rapid timeframe of divergence. However, the earliest components of evolutionary divergence in rapid radiations are often difficult to resolve (Prum et al. 2016;Suh 2016;Moreira and Schrago 2018), and although saturation may play less of a role in younger radiations, other effects may produce similar challenges for their resolution.
As highlighted above, ILS is a major potential source of gene-tree conflict. The ancestral population that gave rise to the mbuna potentially retained substantial genetic polymorphism through multiple early speciation events (Moran and Kornfield 1993). Retained ancestral polymorphisms due to ILS can contribute to substantial gene-tree versus species-tree incongruence in radiations, potentially obscuring or even overriding any molecular phylogenetic signal that arose during their diversification. Methods to efficiently reconstruct species trees in the face of ILS are now well established, making use of coalescent theory by modeling the multispecies coalescent (MSC) process (e.g., Maddison 1997;Rosenberg 2003;Knowles and Carstens 2017;Degnan and Rosenberg 2009;Liu et al. 2010;Bravo et al. 2019;Flouri et al. 2019). At some point, ILS could become insurmountable and render robust phylogenetic reconstructions unachievable, but incorporating ILS into how we model and reconstruct the early branches of explosive radiations like the mbuna will likely help to resolve their phylogeny.
Introgressive hybridization presents another source of gene-tree versus species-tree conflict that can effectively blur phylogenomic signal (Streelman et al. 2004;Mims et al. 2010;Brawand et al. 2014). Hybridization events are departures from strict bifurcation of the species tree, and introduce greater incongruence among gene-trees through the lateral transfer of genes among distinct evolutionary lineages (Braun et al. 2019). Although strict MSC methods account for ILS, they do assume limited introgression among species. Only recently have methods been developed that extend MSC methods to allow for introgression when reconstructing relationships (Degnan 2018;Flouri et al. 2019), but these might not be able to cope with the level of rampant introgression that is assumed to be present in many 2022 SCHERZ ET AL.-PHYLOGENOMICS OF MBUNA CICHLIDS 743 cichlid radiations (Irisarri et al. 2018;Malinsky et al. 2018;Salzburger 2018;Svardal et al. 2021).
In the mbuna, introgression is thought to be exceedingly common, because they can often be readily hybridized in the laboratory and have been documented to hybridize in the Lake (Table 1). Opportunities for such introgressive hybridization are plentiful for mbuna cichlids. Many closely related congeneric lineages often coexist in microsympatry (syntopy), in some cases not just locally, but across the entire distribution of the lineages. This is seen frequently in the morphologically unusual genus Labeotropheus, which has a uniquely hypertrophied snout and highly derived rectangular jaws (Albertson and Pauers 2019;Conith et al. 2019). At numerous sites around Lake Malaŵi, pairs of deep-bodied and shallow-bodied forms occur that are microendemic to single locations (Ribbink et al. 1983a,b;Konings 2007). In each case, the two forms tend to differ not only in body depth and ecology, but also in male nuptial coloration (Ribbink et al. 1983a,b;Konings 2007;Pauers 2010). For instance, at Thumbi West Island in southern Lake Malaŵi, the shallow-bodied Labeotropheus trewavasae and deep-bodied Labeotropheus artatorostris (Pauers 2017) occur in syntopy. As these two species are known to produce F1 hybrids in the lab, and have been implicated in hybridization events at Thumbi West Island, they might be particularly difficult to disentangle genetically. Hence, this species pair offers an exemplary system of syntopic lineages in which to test whether whole-genome data can readily differentiate such sympatric, closely related mbuna species.
In addition to elevated ILS and introgressive hybridization, many recent explosive radiations exhibit little genetic differentiation over the majority of the genome, and contain relatively few phylogenetically informative loci (Campbell and Bernatchez 2004;Via and West 2008;Nosil et al. 2009;Brawand et al. 2014;Kautt et al. 2016;Cai et al. 2021). For example, variation in nuclear genes, which tend to be conserved components of the genome, may be relatively uninformative for phylogenetic reconstruction Brawand et al. 2014). As a result, only a comparatively small portion of the genome can be expected to hold phylogenetically informative loci. Methods based on sequencing only small fragments of the genome (e.g., RADseq, UCEs) are unlikely to yield many informative loci. Whole-genome resequencing, on the other hand, has the potential to make available all of the phylogenetically informative portions of the genome for analysis.
No single method currently exists to identify and overcome the heterogeneous sources of gene-tree conflict and lack of phylogenetic resolution in rapid radiations. As a result, the emerging approach is to integrate and compare results from multiple methodological approaches, each of which targets specific sources of conflict and disparity (Morales-Briones et al. 2021), often interrogating the data in a node-by-node approach (Singhal et al. 2021). So far, most such approaches have worked with reduced representation sequencing data sets, and not the whole genome of the organisms in question. In this study, we adopt an integrative approach to investigate the power of whole-genome resequencing data to resolve the phylogeny of the exceptionally rapid radiation of Lake Malaŵi rock-dwelling cichlids, taking especially ILS and introgression into account. We demonstrate that the resolving power of even wholegenome level data is limited at basal nodes within such an extremely rapid radiation, but, encouragingly, can provide excellent resolution among genera and even differentiate very closely related syntopic Malaŵi species.

Sampling, Genome Resequencing, and Processing
A total of 26 resequenced genomes representing 18 species were included for phylogenomic study (Table S2 of the Supplementary material available on Zenodo at https://doi.org/10.5281/zenodo.5164151). New wholegenome sequences were obtained from 19 adult fish in breeding coloration collected with scuba and barrier nets from Lake Malaŵi in 2010 using permits granted from the Malaŵi Parks Department. High-molecularweight DNA was extracted from fin or muscle tissue from all individuals using commercial kits (QiaGen Dneasy Blood & Tissue Kit) and included an RNase A treatment step. DNA integrity was manually inspected on agarose gels and concentrations were determined on a QuBit fluorometer. Genomic libraries were prepared using Illumina TruSeq DNA Nano kits (Illumina Inc., San Diego, CA, USA) aiming for 350-bp insert sizes. Genomic libraries were then paired-end sequenced (2×150 bp) on a HiSeq 4000 or HiSeq X-Ten Illumina platform at the Beijing Genomics Institute (BGI, Shenzhen, China). Pooling four to five individuals per lane resulted in an average effective genome coverage (counting only reads with mapping quality ≥30, nucleotides with base quality ≥20, and no read duplicates) of approximately 20× per individual.
After demultiplexing, we converted raw reads to unmapped BAM files for long-term storage using Picard tools v.2.9.4 (https://broadinstitute.github.io/picard), adding read group information, and marking adaptor sequences in the process (using the FastqToSamMark and Illumina-Adapters modules). Reads were then converted back into fastq format ( Considering all samples together, we jointly called variants and individual genotypes with freebayes v. 1.1.0 (Garrison and Marth 2012) using default parameters and applying standard quality filters (mapping quality ≥30, base quality ≥20). Subsequently, hard filters were applied using the vcffilter script from the vcflib package (https: to remove low-quality variant sites. Variant representation was normalized using vt norm (Tan et al. 2015). Individual genotype calls based on a read depth smaller than five were set to missing for all downstream analyses. We excluded small scaffolds and restricted analyses to data from the 22 chromosomes only.
After calling single nucleotide polymorphisms (SNPs), VCF files were filtered for the following parameters using vcftools 0.1.15 (Danecek et al. 2011): indels removed, 2 alleles allowed per site, maximally 10% missing data, minimum quality 30, minimum depth 10, and maximum depth 50. After subsequent filtration steps, only variant sites were retained using the -non-ref-ac-any 1 command in vcftools 0.1.15. This filtering resulted in a master VCF file containing 7,421,398 unpruned, genome-wide SNPs. This VCF file was then parsed into separate data sets containing SNPs from noncoding (7,121,933 SNPs) and coding (299,465 SNPs) regions based on CDS annotations in the M. zebra reference genome (Conte et al. 2019).

Data Sets
Genomic data sets comprised of two different taxon sampling schemes were analyzed, hereafter referred to as our Species data set and Population data set: These datasets were deposited on Dryad (https://doi.org/10.5061/dryad.ffbg79cvr).

Species Data Set
For some analyses of our Species data set, SNPs were filtered based on linkage disequilibrium (LD) with bcftools 1.11-88 (Li et al. 2009;Danecek et al. 2021). The +prune plugin was used to calculate r 2 values across sites in windows of 500 kb and exclude SNPs in high LD (r 2 > 0.9). This reduced the size of the noncoding and coding data sets to 495,087 and 38,104 SNPs, respectively. To maximize signal and reduce potential bias due to selection on coding regions, we focused mostly on the noncoding LD-pruned data set for our analyses.
For some analyses, data sets were phased with Beagle 5.1 (Browning et al. 2018), which implements the method of Browning and Browning (2007). We did not specify a genetic panel, as Beagle does not require one to phase accurately, because it generates its own reference panels during analyses (Browning et al. 2021). Beagle assumed a constant recombination rate of 1 cM per 1 Mb to phase genotypes. VCF files were also converted to nexus and phylip files as needed using the Python script vcf2phylip (Ortiz 2019).
To assess genome-wide LD, we plotted LD decay directly from the master VCF file containing both noncoding and coding SNPs using the program PopLDdecay (Zhang et al. 2019) and observed that the decay of LD (r 2 ) leveled off to a mean of ∼0.17 at a distance of roughly 25 kb between sites (Fig. S1 of the Supplementary material available on Zenodo). This distance was relevant in the selection of a sliding window size to produce input "gene" trees for our ASTRAL-III analysis (see below).
Phylogenomics.-We compared phylogenomic trees reconstructed using four methods, including one maximum-likelihood-based concatenated approach, and three coalescence-based methods. We considered the consistency of phylogenetic relationships (reproducibility) across independent methods as an indication of phylogenomic hypothesis robustness (Suh 2016). The non-Lake Malaŵi cichlids H. tweddlei and A. bloyeti were set as outgroups (Malinsky et al. 2018).
Coalescence species tree inferences were performed using three different approaches: 1. SVDquartets (Chifman and Kubatko 2014) implemented in PAUP* 4.0a (Linux build 166) (Swofford 2002). For this analysis, the noncoding LDpruned data set was analyzed using exhaustive quartet sampling and 100 bootstrap replicates. The inferred tree was calculated with the QFM algorithm (Reaz et al. 2014) that optimizes the tree based on the maximum quartet consistency.
2. SNAPP (Bryant et al. 2012) implemented in BEAST2 2.6.3 (Bouckaert et al. 2014). As the analysis of large phylogenomic data sets remains computationally prohibitive with SNAPP, we ran analyses across multiple (n = 10) SNP data sets that each contained a different sampling of 2200 SNPs (n = 100 SNPs randomly sampled from each chromosome of the noncoding LD-pruned data set using bcftools +prune). XML files were generated with the Ruby script snapp_prep.rb (available at https://github.com/mmatschiner/snapp_prep). A starting tree was used as a guide and the ingroup crown age was constrained using a log normal distribution with a mean of 0.5993 and standard deviation of 0.18 based on divergence dating results from Schedel et al. (2019). SNAPP analyses were run in parallel for 5 million generations, sampling every 5000 steps. Log files were examined in Tracer 1.7.1 (Rambaut et al. 2018) and ESS values checked for stationarity across all metrics. After combining all trees sampled from the 10 independent runs in LogCombiner 2.6.3 (discarding a burn-in of 10% from each), a maximum clade credibility tree was produced in TreeAnnotator 2.6.3 assigning median heights at nodes. Maximum clade credibility trees from each of the 10 different runs were also obtained and compared with one another to assess topological convergence between the subsampled data sets.
3. ASTRAL-III 5.6.3 (Zhang et al. 2019) (hereafter "ASTRAL"). For this analysis, we evaluated phased SNPs from the master VCF. Unlike the other three phylogenetic analyses that evaluated the LD-pruned noncoding data set, we here assessed phased SNPs from the genome-wide unpruned VCF file (containing 7,421,398 SNPs). This was done because ASTRAL requires input gene-trees and thus essentially relies on linked data. Maximum likelihood "gene" trees were first generated by conducting RAxML v8.2.12 (Stamatakis 2014) analyses with the GTRCAT model on SNPs that were partitioned into nonoverlapping coordinate windows of 25 kb (-windType coordinate -w 25000) using the Python script raxml_sliding_windows.py (available at https://github.com/simonhmartin/genomics_ general). This window size was selected based on the leveling off of the LD decay (r 2 ) that we observed at approximately this distance (see Fig. S1 of the Supplementary material available on Zenodo). To minimize the generation of suboptimal input trees, maximum likelihood trees were calculated only for windows bearing at least 50 SNPs (-M 50). In total, 26,581 "gene" trees were estimated using this approach. An ASTRAL analysis was then carried out on these trees under default parameters. We also tested for polytomies at deep nodes in the tree with the method of Sayyari and Mirarab (2018), which is based on a Chi-Square test among quartet frequencies for nodes, implemented with the -t 10 command in ASTRAL using these 26,581 trees.
Distance matrices were produced using phased SNPs from the master VCF file and the noncoding and coding LD-pruned data sets with the Python script distMat.py (available at: https://github.com/simonhmartin/geno mics_general) (Table S3 of the Supplementary material  available on Zenodo). Phylogenetic neighbor networks were then constructed from these three matrices with NeighbourNet in SplitsTree 4.16.1 (Hudson and Bryant 2006) to visualize the general structure of our SNP data (Fig. S2 of the Supplementary material available on Zenodo). We also assessed structure among our ingroup species with genomic principal components analysis (PCA). This was done in plink 1.90 (Chang et al. 2015) using the noncoding LD-pruned data set with the two outgroup taxa removed.
Topology weighting. Based on the robustly supported clades recovered using the tree inference methods outlined above (see Results), we conducted fiveclade (outgroup included) and four-clade (outgroup excluded) TWISST (Martin and van Belleghem 2017) analyses to examine whether the basal polytomy within the mbuna was due to regionally conflicting signal across the genome. These clades were those recovered robustly in likelihood and coalescent tree building above. We followed the pipeline of Martin and van Belleghem (2017) using the scripts parseVCF.py, raxml_sliding_windows.py, twisst.py, and plot_twisst.R (available at: https://github.com/simonhmartin/geno mics_general/). Phased SNPs from the master VCF file were analyzed (7,421,398 SNPs for the five-clade analysis and 3,546,543 SNPs for the four-clade analysis excluding the outgroups). A window size of 250 SNPs and the GTRCAT model implemented in RAxML (Stamatakis 2014) were used. Exploratory analyses using 50-and 100-SNP windows showed very little difference among results (data not shown). Haplochromis tweddlei and A. bloyeti were assigned as the outgroup clade in the fiveclade TWISST. In total, 29,675 250-SNP windows were analyzed in the five-clade TWISST and 14,175 in the fourclade TWISST. Results were visualized in R 4.0.2 (R Core Team 2014) using modifications of the scripts provided alongside the TWISST github page.
Analysis of introgression. We calculated Patterson's Dstatistics to assess for introgression across the 16 ingroup mbuna using the ABBA-BABA test as implemented in the program Dsuite ). This test is applied to biallelic SNPs across four taxa and assumes a pectinate tree topology typically given as ([{P1,P2},P3],O). The outgroup (O) helps to define the ancestral allele (A) from the derived allele (B) and site patterns (BBAA, ABBA, and BABA) are counted. Under the null model where only ILS is present (i.e., no gene flow, D = 0), ABBA-BABA patterns are expected to occur in equal frequency, but a significant divergence from this suggests introgression between P3 and either P1 or P2 ). Using the 7,421,398 SNPs from the genome-wide unpruned VCF file, our ASTRAL tree as a guide tree, and A. bloyeti set as the outgroup, we assessed all possible three taxon combinations (560 in total) with the Dtrios function. Each trio was ordered 746 SYSTEMATIC BIOLOGY VOL. 71 according to the tree. Standard jackknife blocks (×20) were used to determine if the resulting D-statistic values differed significantly from zero. To account for multiple tests, P-values were adjusted in RStudio 4.0.3 by applying the false discovery rate method of Benjamini and Hochberg (1995) with the stats package (command: p.adjust(p_values, method = "fdr")). An of 0.01 was applied to conservatively identify statistically significant D-statistic values. To visualize species pairwise comparisons of D-statistic scores, a heatmap was plotted using the Ruby script plot_d.rb (available at: https://github.com/mmatschiner).
After inferring numerous cases of hybridization between the four main clades (A, B, C, and D), we used the Dinvestigate function in Dsuite  to locate specific genomic regions that show elevated patterns of introgression. To examine introgression across the four main mbuna clades, Dstatistic and related metrics (i.e., f d ,f dM , and df) were calculated. These were estimated across sliding windows of 50 informative SNPs at step intervals of 25 SNPs (the program defaults) for trios (240 in total) that yielded significant genome-wide D-statistic values in the Dtrios analysis. The statistics were specifically calculated for trios in which P1 and P2 (i.e., of the ABBA-BABA tests) were assigned to taxa of the same main clade. Astatotilapia bloyeti was used as the outgroup (O). We then evaluated the resulting f dM scores (described in Malinsky et al. 2015), a modified version of the f d statistic that was devised by Martin et al. (2015), to detect introgressed loci. Values of f dM are distributed around zero (no introgression), range from −1 to 1, and quantify the proportion of introgression between P3 and P2 (positive values) and also P3 and P1 (negative values) . We inferred genomic regions to contain introgressed loci if they exhibited f dM scores greater than 0.1 or less than −0.1. These arbitrary thresholds were selected after viewing f dM scores plotted along all chromosomes for several of the trios showing the highest levels of introgression and observing that most values fall within this range and with many obvious narrow peaks extending beyond these thresholds (see Fig. S3 of the Supplementary material available on Zenodo). Taken across all 240 trios, the putative introgressed regions based on these thresholds spanned 21.695% of the 22 chromosomes. The genomic coordinates of these regions were concatenated into a single sorted and merged BED file with bedtools 2.26.0 (Quinlan and Hall 2010), which was subsequently called with bcftools to remove introgressed SNPs from the genome-wide unpruned master and noncoding LDpruned VCF files. This resulted in data sets that included 6,163,040 and 110,402 SNPs, respectively.
To assess the impact of removal of the putative introgressed SNPs on phylogenetic reconstruction and resolution among clades, SVDquartets and IQ-TREE analyses were conducted once more on the "nonintrogressed" noncoding LD-pruned data set with the same parameters as described above in the Phylogenomics subsection. Given the poor performance of SNAPP in our main analysis and the computational constraints it imposes, this program was not used to analyze a "nonintrogressed" data set. An ASTRAL analysis was also performed using the "nonintrogressed" version of the genome-wide unpruned VCF file following the same workflow as outlined before. 21,844 "gene" trees were estimated with RAxML and then analyzed in ASTRAL.

Population Data Set
The population data set was separated from the master VCF file by removing irrelevant individuals, and then removing all invariant sites (-non-ref-ac-any 1) in VCFtools. The resulting data set consisted of 3,731,859 variant SNPs, of which 1,349,540 SNPs were variable within Labeotropheus individuals alone.
Gene flow analyses and topology weighting. We analyzed gene flow among syntopic Labeotropheus species using ADMIXTURE 1.3.0 (Alexander and Lange 2011), F ST , and genomic PCA. ADMIXTURE was performed with k=1-5 groups, and the optimal grouping assessed by comparing the cross-validation (CV) error. Weir and Cockerham F ST was estimated with 100 kb windows in vcftools 0.1.15 (Danecek et al. 2021), with negative F ST values curtailed to 0. Genomic PCA was calculated using plink 1.90b6.16 (Chang et al. 2015), retaining 38,314 variants for the 10 individuals.
To more closely assess the possibility that gene flow has occurred between Labeotropheus and M. zebra, as has been suggested in the past ), we calculated D-statistics in Dsuite , and performed a four-clade TWISST (Martin and van Belleghem 2017) analysis as outlined above, with the clades being L. artatorostris, L. trewavasae, M. zebra, as well as A. bloyeti included as the outgroup, using 250-SNP windows, not thinned for linkage. This window size was selected for consistency with the Species data set TWISST analysis. In total, 14,915 250-SNP windows were analyzed.

Improved but Still Limited Phylogenomic Resolution at the
Base of the Mbuna Cichlids Using SVDquartets on the LD-pruned data sets, we recovered a more robustly supported tree when analyzing the 495,087 SNP noncoding as compared with the 38,104 SNP coding data set (Fig. S4 of the Supplementary material available on Zenodo). Therefore, we focused further analyses on the noncoding SNPs with IQ-TREE and SNAPP and present those results below. The ASTRAL analysis was conducted using "gene" trees generated from the full set of SNPs (coding + noncoding) of the master VCF file. Phylogenomic analyses consistently recovered four main clades within the mbuna: (A) Labeotropheus, (B) Petrotilapia + Tropheops, (C) Melanochromis + Labidochromis + Iodotropheus, and (D) Genyochromis + Cynotilapia + Maylandia + Chindongo (Fig. 1). Genomic PCA also separated these clades into four distinct clusters (Fig. S5 of the Supplementary material available on Zenodo). All reconstructions agree in placing Iodotropheus within Labidochromis with full support. They also agree in placing Melanochromis as the sister of these taxa with full support. Ignoring the position of the root, SVDquartets and IQ-TREE recovered the same quartet topology among the four clades of A,D|B,C, whereas SNAPP and ASTRAL recovered the quartet topology A,C|B,D, but with different placements of the root (Fig. 1). The third possible quartet, A,B|C,D, was not recovered by these analyses. SVDquartets and IQ-TREE trees yielded almost identical overall topology, differing only in the monophyly of Tropheops species in Clade B (monophyletic in IQ-TREE, paraphyletic in SVDquartets), and in the branching order of Genyochromis, Cynotilapia, and Chindongo in Clade D. The latter three taxa branched in that order in all analyses except IQ-TREE, where Chindongo fell sister to Maylandia, and the other two (Genyochromis + Cynotilapia) as a sister clade (Fig. 1).
Labeotropheus (Clade A) was found to be sister to all other clades in ASTRAL, SVDquartets, and IQ-TREE, whereas Clade C was recovered as the earliest diverging mbuna lineage in the SNAPP analysis. Three of the 10 SNAPP runs conducted on the independent sets of 2200 SNPs recovered Labeotropheus at the base of the radiation (data not shown). Clade C was recovered as monophyletic in only five of these independent SNAPP analyses. The maximum clade credibility tree found by combining these analyses shows Clade C as monophyletic albeit with poor support. We think the poor resolution captured across the SNAPP analyses may be due to the limited number of SNPs used per run. Despite high ESS values, the 10 analyses failed to individually converge on a common topology, and we therefore view the maximum clade credibility tree obtained from evaluating trees taken across all 10 analyses as potentially spurious.

Sources of Poor Resolution at the Base of the Mbuna Cichlids
We identified several factors that likely contributed to the poor resolution of the nodes at the base of the mbuna tree. First, despite moderately high support, a polytomy could not be rejected for one of the deepest nodes within the mbuna (that at the split of Clades B + D: 2 test with the method of Sayyari and Mirarab (2018), P = 0.2048) in the ASTRAL tree (Fig. 1). All other nodes across that tree remained robust, rejecting the null hypothesis of polytomy at = 0.05. Second, in our IQ-TREE analysis, site concordance factors were low at several of the basal nodes, although SH-aLRT tests produced values of 100 across all nodes (Fig. 1).
Examining quartet frequencies at the nodes of our ASTRAL tree, we were surprised to find that there were no nodes at which the assumption of no ILS was obviously violated. In all cases, the second and third hypotheses had similar frequencies (under introgression, we would expect unequal frequencies of these alternative hypotheses). To investigate this further, we conducted an analysis of D-statistics across all taxa using Dsuite  on SNPs from the master VCF file. Significant D-statistic values (P < 0.01) potentially indicative of low levels of introgression were inferred in 422 of the 560 possible trios and between 96 unique mbuna species pairs (Fig. 2a, Table S4 of the Supplementary material available on Zenodo). Tropheops tropheops in particular exhibited many patterns of introgression with taxa from Clades A, C, and D. In line with hybridization that has been previously reported (Table 1), M. zebra and the two Labeotropheus species also show signal of potential introgression with other mbuna.
A neighbor network constructed from the genomic distances among taxa (Fig. 2b, Fig. S2 and Table S3 of the Supplementary material available on Zenodo) also recovers the four clades outlined above, as well as many of the subclades. The position of Iodotropheus within Labidochromis, the close relationships of the other genera, and the very deep divide between the Tropheops species are all notable and help to explain the inconsistency of their relationships among reconstruction methods. The central knot of this network contains little structure, as already indicated by the lack of consistency and support in these nodes in ML and coalescence reconstructions, and inability to reject a basal polytomy in ASTRAL analyses.
Given the broad extent of hybridization inferred here across the mbuna, we located genomic regions exhibiting elevated patterns of introgression and masked SNPs from these loci in follow-up phylogenomic analyses to gauge their impact on clade resolution. The putative introgressed SNPs that were identified across analysis of 240 trios spanned 21.695% of the 22 chromosomes (165,098,567 of 760,997,612 bp). After removal of these SNPs, the unrooted quartet branching pattern of the main clades persisted in both the SVDquartets and ASTRAL analyses (A,D|B,C vs. A,C|B,D, respectively), but was inconsistent between the IQ-TREE analyses (Fig. 1, Fig. S6). With IQ-TREE, we recovered the quartet topology A,C|B,D from the introgressed-SNPpruned data set despite observing A,D|B,C in the original analysis. Although the splitting pattern among the main clades remained unresolved after excluding putative introgressed SNPs, support for the monophyly of Clade D drastically improved in subsequent IQ-TREE, SVDquartet, and ASTRAL analyses. Support for Clade B, on the other hand, was reduced in each of these additional analyses. Having recovered four mbuna clades consistently in MSC and concatenated approaches as well as in the distance-based neighbor network, we were interested to know whether specific regions of the mbuna genomes were giving conflicting signals with respect to the deep branches. We investigated this using topology weighting analysis (TWISST) based on these clades and ran the analysis with and without outgroups. Our analysis without outgroups (containing 3,546,543 SNPs) found nearly identical weighting between the two quartets A,C|B,D (the topology found by ASTRAL and SNAPP) and A,B|C,D (the topology not recovered by any of our reconstructions). Both of these were weighted stronger than A,D|B,C (the topology recovered by SVDquartets and IQ-TREE; Fig. 3b; Fig. S7 of the Supplementary material available on Zenodo). There was little indication of particular genomic regionalization of quartet weighting (as exemplified in Fig. 3c). The analysis with outgroups recovered the highest average weighting for the topology R((A,B),(C,D)), slightly over R((A,C),(B,D)) (root here indicated with R) (Fig. 3d; Fig. S8 of the Supplementary material available on Zenodo). These two topologies had substantially higher average weighting than the next best topology (Fig. 3e), but none of the three highest weighted topologies correspond to those recovered in our phylogenomic reconstructions. Weighting was extremely heterogeneous across the genome based on 250-SNP windows, with two adjacent windows seldom favoring the same topology (Fig. 3f). Together, the results of our TWISST analyses show highly heterogeneous signal across the genomes of these fishes, which corresponds to the lack of clear resolution found in other analyses.

Genomic Distinction of Two Syntopic Labeotropheus Species
We found robust support for the genomic distinction and reciprocal monophyly of the closely related and syntopically occurring L. artatorostris and L. trewavasae. The two species have a genome-wide weighted Weir and Cockerham F ST of 0.2353 (mean = 0.1226; Fig. 4a). In a genomic PCA, the two species were strongly separated in PC1, explaining 35.4% of the genomic variance (Fig. 4c). The remaining PCs each separated a single individual from the rest, as seen in Fig. 4c for PC2 (other PCs not shown). ADMIXTURE with K = 2 divided the two species, although K = 1 had the lowest CV error, possibly due to our small sample sizes per population (Fig. 4d,e).
We also found little evidence for the introgression between Labeotropheus species and M. zebra that has previously been hypothesized . ABBA-BABA tests yielded a D-statistic of 0.031 (P = 5.251× 10 −14 ; BBAA 109,799, ABBA 25,628, BABA 24,090.4), f4-ratio 0.0101. TWISST showed dominant signal for monophyly of Labeotropheus, with very little indication of introgression from M. zebra into the genome, though marginally more into L. trewavasae than into L. artatorostris ( Fig. 4f-h; Fig. S9 of the Supplementary material available on Zenodo).

DISCUSSION
The integrative and stepwise approach taken here, based on SNP data sets from complete genome sequences of 16 species in 10 genera, was able to provide clear resolution for some, but not all, of the phylogenetic relationships among the Malaŵi mbuna cichlids. With the most comprehensive phylogenomic sample of mbuna cichlids to date, we were able to recover four main clades within the mbuna with high support. Most previous works have included exemplars of each of these clades but lacked the data to test their interrelationships (Table S1 of ]), these are likely to fall within these four main clades. As phylogenomic data sets frequently yield spuriously high statistical support (Rokas and Carroll 2006;Singhal et al. 2021), the consistent recovery across reconstruction methods (Fig. 1) provides greater confidence than any single method alone in the existence of the four main clades.  However, the resolution among these four clades at the base of the mbuna phylogeny remains uncertain. Our phylogenomic analyses provided moderate to high support for two alternative quartet configurations among the clades, and mostly agreed in placing Clade A, consisting only of Labeotropheus, sister to all other sampled mbuna (Fig. 1, Fig. S6 of the Supplementary material available on Zenodo). However, the branches separating Clades B, C, and D were very short, received low support, and failed the polytomy test implemented in ASTRAL. Paradoxically, the one quartet configuration that was not recovered by any of our coalescent or maximum likelihood reconstruction methods, A,B|C,D, was found to be the most frequent local topology 752 SYSTEMATIC BIOLOGY VOL. 71 across the genome in our topology weighting analysis (Fig. 3). Hybridization and ILS both likely limit our ability to resolve the evolutionary history of the mbuna to some degree. We were unable to unambiguously resolve the interclade relationships even after excluding putative introgressed SNPs in subsequent phylogenomic analyses. Following the pruning of introgressed loci, SVDquartets, ASTRAL, and IQ-TREE analyses still failed to converge on a common topology (Fig. S6 of the Supplementary material available on Zenodo). Both ASTRAL and IQ-TREE recovered a quartet pattern of A,C|B,D, but SVDquartets reported A,D|B,C, albeit with weak support. Although resolution of the backbone is slightly improved when putatively introgressed loci are excluded with SVDquartets and ASTRAL analyses, support across these deep nodes is diminished in the IQ-TREE phylogeny. Furthermore, the internodes in our ASTRAL tree separating the main four clades remained very short after controlling for introgression. Taking all this into account, it seems that the basal nodes within the mbuna radiation defy resolution even with methods accounting for ILS and introgression from whole-genome data.
On the other hand, despite the polytomy at the base of the radiation and the short timeframe over which the mbuna are thought to have diversified, the entire mbuna phylogeny is not intractable. Within each of the four robust clades, we were able to achieve good phylogenomic resolution, with only little inconsistency among reconstruction methods for most relationships. Although the monophyly of Tropheops and Labidochromis are called into question (see below), there was substantial phylogenetic structure recovered among the disparate lineages we sampled ( Figs. 1 and 2), and our sampling of couplets of species within six of the genera generally provided unambiguous monophyletic relationships. This robust ability to resolve these relationships runs contrary to the expectation that congeneric Malaŵi cichlids would be difficult to distinguish genotypically (Moran and Kornfield 1993;Albertson et al. 1999;Hulsey et al. 2010).
This ability of the whole-genome SNP data to resolve the more recently diverged mbuna relationships was particularly reinforced by the reciprocal monophyly, elevated F ST , and lack of extensive admixture that we found in the two sympatrically occurring Labeotropheus species (Fig. 4), despite the fact that they hybridize in the lab (Table 1) and can often be observed to feed side-by-side on the same rock (Hulsey pers. obs.). The data showed that these species have only low levels of introgression, but that other species, especially Tropheops and taxa belonging to Clade D, have much higher levels across the mbuna radiation (Fig. 2a). Our results suggest that interclade hybridization is evidently common within the mbuna cichlids, but may not pose a major threat to the reconstruction of a largely bifurcating phylogeny of these fishes, in contrast to previous assumptions Malinsky et al. 2018;Svardal et al. 2021). Although the removal of introgressed loci failed to resolve the basal branching order of the four main clades with certainty, it improved resolution among other parts of the tree, particularly for the monophyly of Clade D and relationships within it. This action however, resulted in a significant reduction of data set size (e.g.,110,402 SNPs vs. 495,087 SNPs for the noncoding LD-pruned data set) and as a consequence, the support for some relationships was diminished. As a whole, these results suggest that genome-wide SNP data in the mbuna and likely other recent radiations might commonly contain enough informative markers to account for introgression and still recover underlying phylogenetic history.
Due to their intermediate morphology, it has been suggested that Tropheops might be the result of hybrid speciation between Labeotropheus and Maylandia ). This hypothesis is rejected by our phylogenomic reconstructions, which robustly place Tropheops in a clade with Petrotilapia, whereas Maylandia is more closely related to Genyochromis, Cynotilapia, and Chindongo. However, we did find evidence for considerable gene flow between Tropheops and Labeotropheus, Chindongo, and Genyochromis (Fig. 2a). Although this does not directly support a hybrid origin of Tropheops, it certainly does warrant further study. SVDquartets recovered Tropheops as paraphyletic with respect to Petrotilapia, but no other method supported this arrangement. Our phylogenomic reconstructions also reveal that Iodotropheus falls within Labidochromis. This is strong evidence that the two genera are synonymous (I. sprengerae is the type species of Iodotropheus). However, formal taxonomic synonymization of these genera should be forestalled until a more comprehensive study is undertaken on the two other valid Iodotropheus species and 16 other valid Labidochromis species, including the type species of Labidochromis, L. vellicans (Fricke et al. 2020). Curiously, Malinsky et al. (2018) did find Maylandia to be closely associated with Petrotilapia and Tropheops, differing substantially from our findings. However, their M. zebra specimen was a different individual than ours, so it may be that the two sampled individuals assigned to this species were not conspecific. Also, Labeotropheus was missing from most of their analyses. In our case, the incorporation of a second Maylandia species, M. xanstomachus, helps to inform and support the relationships that we recovered. Additional intergeneric sampling could further clarify and reinforce our general understanding of mbuna relationships.
Resolution of the base of the mbuna radiation was often one of the primary foci of earlier phylogenetic studies (Moran and Kornfield 1993;Albertson et al. 1999;Hulsey et al. 2017). However, based on our large genome-wide data set, it seems likely that the lack of resolution reflects biological processes that render this part of the mbuna phylogeny highly intractable to systematic reconstruction. In contrast, reconstruction of more recent bifurcation events in the mbuna radiation were readily resolved, and future phylogenetic studies 753 should focus primarily on these more recent, tractable, and diagnosable evolutionary relationships. The reconstruction of a phylogeny and subsequent species delimitation of mbuna genera across the whole of Lake Malaŵi may therefore be possible with whole-genome sequences. Resolving the age-old issue of how many species of cichlids occur in one of the most diverse radiations of fishes in the world may now be feasible.

Between a Rock and a Hard Polytomy
The major clades of mbuna recovered may represent the product of a phylogenetically unresolvable burst of diversification. This type of rapid divergence is thought to characterize many adaptive radiations (Suh 2016;Moreira and Schrago 2018;Pease et al. 2018;Braun et al. 2019;Cai et al. 2021;Gagnon et al. 2021;Hime et al. 2021;Morales-Briones et al. 2021;Singhal et al. 2021). The power to resolve some polytomies may be proportional to the volume of relevant data available (i.e., more loci should better resolve such recalcitrant nodes). In the past, the nodes lacking resolution at the base of the mbuna were treated as soft polytomies that could be resolved with sufficient data (e.g., Hulsey et al. 2017). However, our results show that sufficiently rapid radiations can simply defy phylogenomic resolution, even with whole-genome resequenced data. The relationships among the major clades of mbuna could not be confidently resolved, and the relationships agreed only between two of our four phylogenomic reconstructions (IQ-TREE and SVDquartets) obtained through the main analyses (Fig. 1). Although we would welcome innovative attempts to resolve these relationships, our data call into question whether these early mbuna relationships will ever be resolved with molecular data; this may be a hard polytomy.
Our ability to resolve rapid diversification events is often dictated by the interactions of evolution with time. In some cases, extinction of lineages involved in a burst of diversification increases the remaining internal branch lengths, potentially making early nodes within a radiation easier to resolve as clades age. However, saturation of molecular data, especially when lineages do not undergo substantial extinction, diminishes and eventually extinguishes phylogenetic signal (Xia et al. 2003). The number of orthologous SNPs also decreases over time especially in noncoding regions of the genome, which reduces the usability of SNP-based data sets to resolve diversification bursts of increasing age (Leaché and Oaks 2017). Gene duplications and loss, as well as genomic rearrangements can also pose problems (Bravo et al. 2019). Together, these processes decrease the power of genomic data to resolve bursts of diversification of increasing age and thereby result in decreasing power to resolve polytomies over time, irrespective of whether they are the result of multifurcation or rapid sequential bifurcation.
However, because the mbuna radiation is extremely young (<2 million years; Schedel et al. 2019;Matschiner et al. 2020), especially compared with some other examples of polytomies that are well established in the literature (e.g., that at the base of the Neoaves, ∼66 million years; Suh 2016; Braun et al. 2019), it probably has not been affected significantly by time-related signal loss. Both the short branches recovered by ASTRAL and IQ-TREE and the low support values associated with these nodes across all of our reconstructions are consistent with a very rapid radiation at the origin of the mbuna clade. This has likely contributed to extremely high levels of shared standing genetic variation among species and ILS, making its basal branching relationships difficult to disentangle . The high levels of standing genetic variation may have increased the level of plasticity in the early radiation, and thereby actually facilitated the diversification of the clade by subsequent "genetic assimilation" (see Meyer 1987;Schneider and Meyer 2017). Similar problems have been identified in other genomic analyses of cichlid relationships (e.g., Browning et al. 2018;Salzburger 2018;Olave and Meyer 2020;Svardal et al. 2021). Additionally, our analysis shows that gene flow has been moderately common among clades, consistent with previous reconstructions of this and other African Great Lake cichlid radiations (Irisarri et al. 2018;Malinsky et al. 2018;Salzburger 2018;Svardal et al. 2021). However, this introgression has not been as obfuscating as we had originally expected. Instead of an overall pattern of high introgression, a few specific lineages appear to be the main drivers of introgressive hybridization (Fig. 2a), and most of these involved interclade, and not intraclade, gene flow. Removing those regions that are particularly strongly introgressed resulted in better support for the placement of some introgressing taxa, but did not improve the interclade relationships. So, although it has likely contributed to uncertainty together with other sources of ILS, introgression does not seem to be the main source of the basal node instability in the mbuna.

CONCLUSION AND OUTLOOK
In the Lake Malaŵi mbuna cichlid radiation, we have found that integrating results from four different phylogenomic reconstruction methods that utilized SNP data from whole-genome resequencing provided strong phylogenomic resolution at most levels, but not the basal-most nodes. This provided the most robust phylogenetic hypothesis for this lineage to date, with four major clades being resolved that had only partially been recovered before, and Labeotropheus emerging as perhaps the sister to all other sampled mbuna. We suggest that the short branches of this rapid radiation may represent an unresolvable polytomy at the base of the mbuna. If this is the case, future work with more extensive taxon sampling will not resolve the basal nodes further (Suh 2016), but molecular examinations of other parts of the tree will likely be more fruitful. This also has implications for the taxonomic efforts to describe mbuna diversity. Taxonomists should work to incorporate molecular data, 754 SYSTEMATIC BIOLOGY VOL. 71 especially large genomic data sets, as these efforts should now be seen as finally capable of providing robust and potentially complementary tests of more traditional morphology-based taxonomic diagnoses (e.g., Hanssens 2004;Li et al. 2016;Oliver 2018).
Nevertheless, as we start to accept hard polytomies as describing the evolution of some major clades, it is becoming more important to incorporate these evolutionary patterns into evolutionary studies such as trait evolution. The interest in the evolutionary links between ecology, phenotypes, and genotypes in rapidly diversifying groups is only increasing. Randomization techniques and integration of results across a large set of bootstrapped trees is common practice in current comparative analyses (Losos 1994;Hulsey et al. 2007), but fails to account for the inferential limitations of the underlying data sets (Smith et al. 2020). Moreover, the study of adaptive divergence in explosive radiations has increasingly shown that the evolutionary history of traits themselves sometimes do not correspond to the species tree (Kratochwil et al. 2018;Kautt et al. 2020), so the importance of incorporating ILS into our understanding of trait evolution is likely only going to increase. The use of gene-trees instead of bootstrap trees might make the iterative modeling of evolution over troublesome nodes more informative (Hahn and Nakhleh 2016), but ultimately, what is needed are methods to translate node uncertainty into the models underlying comparative analyses (Singhal et al. 2021). As we continue to improve the methodological toolkit for inference of accurate evolutionary histories, and identify also the limitations of even the most comprehensive genomic data sets, we will be better suited to extract what can and cannot be learned about organismal evolution within the framework of molecular systematics.