Abstract
Free full text
Widespread introgression across a phylogeny of 155 Drosophila genomes
SUMMARY
Genome-scale sequence data have invigorated the study of hybridization and introgression, particularly in animals. However, outside of a few notable cases, we lack systematic tests for introgression at a larger phylogenetic scale across entire clades. Here we leverage 155 genome assemblies, from 149 species, to generate a fossil-calibrated phylogeny and conduct multilocus tests for introgression across nine monophyletic radiations within the genus Drosophila. Using complementary phylogenomic approaches, we identify widespread introgression across the evolutionary history of Drosophila. Mapping gene-tree discordance onto the phylogeny revealed that both ancient and recent introgression has occurred across most of the nine clades that we examined. Our results provide the first evidence of introgression occurring across the evolutionary history of Drosophila and highlight the need to continue to study the evolutionary consequences of hybridization and introgression in this genus and across the Tree of Life.
In Brief
Suvorov et al. use multi-locus data from 155 Drosophila genome assemblies to reconstruct the phylogeny of the genus and estimate divergence times. They used multiple phylogenomic methods to shed light on patterns of gene flow within the genus. Their results suggest multiple instances of introgression across the Drosophila Tree of Life.
INTRODUCTION
The extent of gene exchange in nature has remained one of the most hotly debated questions in speciation genetics. Genomic data have revealed that introgression is common across taxa, having been identified in major groups such as fungi1–3, vertebrates4–7, insects8–10, and angiosperms11,12. The evolutionary effects of introgression are diverse, and are determined by multiple ecological and genomic factors13,14. Once thought to be strictly deleterious, it has become increasingly clear that introgression can serve as a source of genetic variation used during local adaptation15,16 and adaptive radiation17,18. While our understanding of introgression as a widespread phenomenon has clearly improved, it remains unclear how often it occurs across taxa. Ideally, determining the frequency of introgression across the Tree of Life would leverage the signal from systematic analyses of clade-level genomic data without an a priori selection of taxa known to hybridize in nature.
At the phylogenetic scale, hybridization has typically been explored at relatively recent timescales. For example, studies of hybridization between cats (Felidae; 10–12 My; ~40 species19), butterflies (Heliconius; 10–15 My; 15 species8), cichlid fishes from the African rift lakes (0.5–10 My; ~27 species18,20,21), and wild tomatoes (Solanum; ~4 My; ~20 species12) all rejected a purely bifurcating phylogenetic history. In each of these systems introgression has occurred relatively recently, as the common ancestor for each species group occurred no more than 15 million years ago. However, there are also notable exceptions, and evidence for introgression has been found across much deeper phylogenetic timescales within vascular plants11 and primates7. In some species, there is also evidence that introgression has been a source of adaptive genetic variation that has helped drive adaptation (e.g. refs. 2,22–25). These results show how introgression has both (1) occurred in disparate taxonomic groups and (2) promoted adaptation and diversification in some. Notwithstanding key examples4–7,11,12, we still require systematic tests of introgression that use clade-level genomic data that spans both deep and shallow phylogenetic time to better understand introgression’s generality throughout evolution.
Species from the genus Drosophila remain one of the most powerful genetic systems to study animal evolution. Comparative analyses suggest that introgression might be common during speciation in the genus26. Genome scans of closely related drosophilid species have provided evidence of gene flow and introgression9,10,27–32. There is also evidence of contemporary hybridization33–35 and stable hybrid zones between a handful of species36–38. These examples of hybridization and introgression show that species boundaries can be porous but cannot be taken as prima facie evidence of the commonality of introgression. We still lack a systematic understanding of the relative frequency of hybridization and subsequent introgression across Drosophila. Here we analyze patterns of introgression across a phylogeny generated using 155 whole genomes derived from 149 species of Drosophila, and the genomes of four outgroup species. These Drosophila species span over 50 million years of evolution and include multiple samples from nine major radiations within the genus Drosophila. We used two different phylogenetic approaches to test whether introgression has occurred in each of these nine radiations. We found numerous instances of introgression across the evolutionary history of drosophilid flies, some mapping to early divergences within clades up to 20–25 Mya. Our results provide a taxonomically unbiased estimate of the prevalence of introgression at a macroevolutionary scale. Despite few known observations of current hybridization in nature, introgression appears to be a widespread phenomenon across the phylogeny of Drosophila.
RESULTS
A high-confidence phylogeny of 155 Drosophila genomes
We first used genome-scale sequence data to infer phylogenetic relationships among species in our data set. To achieve this, we annotated and generated multiple sequence alignments for 2,791 Benchmarking Universal Single-Copy Orthologs (BUSCOs; v339,40) across 155 independently assembled Drosophila genomes together with four outgroups (3 additional species from Drosophilidae and Anopheles gambiae). We used these alignments, totalling 8,187,056 nucleotide positions, and fossil calibrations to reconstruct a fossil-calibrated tree of Drosophila evolutionary history. Note that the inclusion of Anopheles as an outgroup allowed us to include a fossil of Grauvogelia, the oldest known dipteran, in our fossil calibration analysis, along with several Drosophilidae fossils and/or geological information (i.e., formation of the Hawaiian Islands; Data S1).
Our phylogenetic analyses (see STAR Methods for details) using both maximum-likelihood (ML using the IQ-TREE package) and gene tree coalescent-based (ASTRAL) approaches with DNA data revealed well-supported relationships among nearly all species within our dataset. Phylogenies inferred using these two approaches only differed in three relationships (Figure S1): (i) D. villosipedis was either recovered as sister species to D. limitata + D. ochracea (ML topology) or as a sister to D. limitata + D. ochracea + D. murphyi + D. sproati (ASTRAL topology); (ii) D. vulcana and D. seguyi form monophyletic lineage sister to the D. nikananu + D. spp. aff. chauvacae + D. burlai + D. bocqueti + D. bakoue clade (ML topology) or have paraphyletic relationships where D. vulcana is sister to the D. nikananu + D. spp. aff. chauvacae + D. burlai + D. bocqueti + D. bakoue clade (ASTRAL topology); (iii) D. simulans was recovered as sister either to D. mauritiana (ML topology) or D. sechellia (ASTRAL topology, the latter of which is perhaps more likely to be the true species tree according to an analysis examining low-recombining regions, which are less prone to ILS41. The nodal supports were consistently high across both ML (Ultrafast bootstrap (UFBoot) = 100, an approximate likelihood ratio test with the nonparametric Shimodaira–Hasegawa correction (SH-aLRT) = 100, a Bayesian-like transformation of aLRT (aBayes) = 1) and ASTRAL (Local posterior probability (LPP) = 1) topologies with the exception of D. limitata + D. ochracea + D. villosipedis (UFBoot = 9, SH-aLRT = 81, aBayes = 1) and D. carrolli + D. rhopaloa + D. kurseongensis (UFBoot = 81.2, SH-aLRT = 81, aBayes = 1) on the ML tree, and D. limitata + D. ochracea + D. murphyi + D. sproati (LPP = 0.97) and D. sulfugaster bilimbata + D. sulfugaster sulfurigaster (LPP = 0.69) on the ASTRAL tree. Thus, the phylogeny we report here is the first of the genus Drosophila with almost all nodes resolved with high confidence—recent estimates of the Drosophila phylogeny lacked strong support throughout all tree depth levels42–44.
Erroneous orthology inference as well as misalignment can impede accurate phylogenetic inference and create artificially long branches45. Repeating our ASTRAL analysis after removing outlier long branches via TreeShrink45 resulted in an identical tree topology with the aforementioned ASTRAL tree (Figure S1). Furthermore, an ML topology estimated from the dataset with more closely related outgroup species (see STAR Methods) results in an identical topology with the aforementioned ML tree (Figure S1). The inferred phylogeny from the protein supermatrix showed only four incongruencies with the phylogeny that was inferred from DNA data (Figure S1): (i) D. villosipedis was recovered as a sister species to D. limitata + D. ochracea + D. murphyi + D. sproati; (ii) D. watanabei + D. punjabiensis is sister to the clade containing D. bakoue and D. jambulina; (iii) D. vulcana and D. seguyi show paraphyletic relationships; (iv) Z. vittiger and Z. lachaisei show sister species relationships. We performed further assessment of nodal support with Quartet Sampling11, using the Quartet Concordance (QC) and Quartet Differential (QD) scores to identify quartet-tree species-tree discordance (STAR Methods). At some nodes, an appreciable fraction of quartets disagreed with our inferred species tree topology (QC < 1), and in most of these cases this discordance was skewed toward one of the two possible alternative topologies (i.e. QD < 1 but > 0) as is consistent with introgression. We formally explore this pattern below.
In order to estimate divergence times across the Drosophila phylogeny, we developed five calibration schemes (A, B, C, D and “Russo”; Data S1) used in MCMCtree46 and one scheme based on the Fossilized Birth-Death (FBD) process47 used in BEAST248 (BEAST2 FBD; Data S1). Overall, four of the five MCMCtree schemes yielded nearly identical age estimates with narrow 95 % credible intervals (CI), whereas scheme “Russo” (a fossil calibration strategy closely matching that from43) showed slightly older estimates (Figure S2) with notably wider 95% CIs. Throughout this manuscript we use the time estimates obtained with scheme A. This calibration analysis estimated that extant members of the genus Drosophila branched off from the other Drosophilidae (Leucophenga, Scaptodrosophila and Chymomyza) ~53 Mya (95% CI: 50 – 56.6 Mya) during the Eocene Epoch of the Paleogene Period (Figure 1). The same analysis inferred that the split between the two major lineages within Drosophila—the subgenera of Sophophora and Drosophila—occurred ~47 Mya (95% CI: 43.9– 49.9 Mya; Figure 1; “A” and “B” clades, respectively); previously published estimates of this time include ~32 Mya (95% CI: 25–40 Mya49), ~63 Mya (95% CI: 39–87 Mya50), and ~56 Mya (95% CI not available43). We also note that our divergence time estimates of the Drosophila subgenus (~34 Mya, 95% CI: 31.6 – 36.8 Mya; clades 6 through 9) are somewhat younger than ~40 Mya, a previous estimate reported in51, although the latter had fairly wide confidence intervals (95% CI: 33.4 – 47.6 Mya). On the other hand, divergence time estimates produced by the FBD scheme in BEAST2 tend to be older especially for deeper nodes (Figure S2). Also, CIs estimated by BEAST2 were wider than those from MCMCTREE. This can be explained by the fewer assumptions about fossil calibration placement and age prior specification for methods that rely on the FBD process. Additionally, we note that not all parameters of the BEAST2 FBD calibration scheme converged (i.e., effective sample size < 100) even after 6 × 108 MCMC generations. Thus, the lack of a thorough fossil record within Drosophila makes it difficult to accurately and precisely estimate divergence times, and point estimates of divergence times should be interpreted with caution.
The blue distributions at each divergence point on the tree represent nodal age posterior probabilities from MCMCTree. Grauvogelia and Oligophryne fossils were used to set priors on the age of the root of the tree, Phytomyzites and Electrophortica succini were used for priors for the root of the Drosophilidae family, and Electrophortica succini and Scaptomyza dominicana were used to set priors for the crown group “Scaptomyza”, i.e. Most Recent Common Ancestor (MRCA) node of the Scaptomyza species (scheme A; Data S1). The numbered red circles denote clades for which analyses of introgression were performed. Inset: the phylogenetic and temporal relationships between our distant outgroup Anopheles gambiae, more closely related outgroup species of Drosophilidae (Leucophenga varia, Scaptodrosophila lebanonensis and Chymomyza costata), and the Drosophila genus. A and B denote the two inferred major groups within Drosophila. See also Figures S2, S5 and Data S3.
Widespread signatures of introgression across the Drosophila phylogeny
To assess the prevalence of introgression across the Drosophila tree, we subdivided species into nine monophyletic lineages (herein referred to as clades 1 through 9; Figure 1) and tested for introgression within each clade. These clades correspond to the deepest divergences within the genus, with most having an MRCA during the Paleogene. Clades 4 and 5 are the two exceptions, splitting from an MRCA later in the Neogene. Within each of the nine clades, the MRCA of all sampled genomes ranged from ~10 Mya (Figure 1; clade 2) to ~32 Mya (Figure 1; clade 1). We note that Hirtodrosophila duncani, Drosophila busckii and Drosophila repletoides were not included in these clade assignments as each of these species was the only sampled descendent of a deep lineage; additional taxon sampling is required to assign them to specific monophyletic species groups that could be tested for introgression.
We tested for introgression within each of these nine clades using two complementary phylogenomic methods that rely on the counts of gene trees inferred from the BUSCO loci that are discordant with the inferred trees (hereafter referred to as the discordant-count test or DCT) and the distribution of branch lengths for discordant gene trees (hereafter termed the branch-length test or BLT), respectively, among rooted triplets of taxa (Figure 2). These methods leverage information contained across a set of gene trees to differentiate patterns of discordance that are consistent with introgression from those that can be explained by incomplete lineage sorting alone (see STAR Methods). We found at least one pair of species with evidence of introgression in 7 of the 9 clades according to both DCT and BLT (i.e. the same pair of species showed evidence for introgression that was significant for both tests in the same triplet at an FDR-corrected P-value threshold of 0.05). In clades 1 and 3 there were no species pairs for which the DCT and BLT were significant in the same triplet and both suggest the same introgressing species pair (Data S2). However, both clades had several pairs that were significant according to one test or the other (Data S2). We found even stronger support for introgression using two existing software methods: QuIBL (Data S2), which examines the branch-length distributions of all three gene tree topologies for a triplet8, and HyDe (Data S2), which tests for introgression by counting quartet site patterns52. Specifically, QuIBL detected introgression in 120 of 152 (78.9%) of species pairs detected by both DCT and BLT, as well as 894 additional species pairs not detected by DCT-BLT; we note that BLT and QuIBL approaches are not fully independent, since they both utilize branch-length information. Similarly, HyDe detected introgression in 142 of 152 (93.4%) of species pairs detected by both DCT and BLT, and 898 additional species pairs (the results of HyDe were not qualitatively affected if a more distantly related outgroup, i.e. Anopheles gambiae, was selected, see Data S2). However, we focus here on the intersection between DCT and BLT methods (after correcting each for multiple testing), as this provides a more conservative estimate of the extent of introgression. Supporting this claim, we applied these tests to a gene tree dataset simulated under high levels of ILS53 and observed low false positive rates: 0.054 for DCT, 0.089 for BLT, and 0.009 for their intersection.
(A) The Discordant Count Test (DCT) identifies introgression where a given triplet within the tree shows an excess of gene trees that support one of the two possible divergent topologies. Note that concordant gene trees and corresponding probabilities are included for completeness, although these are not used by our test. (B) The Branch Length Test (BLT) identifies introgression where branch lengths of gene trees that support introgression are shorter than branch lengths of those that support the species tree and the less frequent divergent topology (i.e., the discordant topology putatively due to ILS).
We carried out several analyses to assess the robustness of our results to data quality and evolutionary rate. Specifically, we tested whether aliment length and quality, karyotype evolution 54, and positive selection may potentially affect introgression inference by the DCT and BLT approaches by repeating our analysis after filtering our data to account for each of these potential confounders (STAR Methods). Although these filtering schemes overall reduce the number of introgression events, our major conclusions about introgression within Drosophila still hold (Table S1). Indeed, most of the reduction in signal that we observe appears to be driven by the loss of power due to the reduced size of the filtered dataset (see effect of filtering random sets of gene trees; Table S1), rather than data artifacts and/or biological confounding factors which appear to have only a modest impact.
The number of species pairs that show evidence of introgression in our initial DCT-BLT analysis is not equivalent to the number of independent introgression events among Drosophila species. This is because gene flow in the distant past can leave evidence of introgression in multiple contemporary species pairs. For example, we found evidence for introgression between D. robusta and all five species within the D. americana-D. montana group (see clade 7 in Figure 3). Rather than five independent instances of introgression between species, this pattern could reflect introgression between ancestral taxa that subsequently diverged into the contemporary species. More generally, cases where multiple introgressing species pairs each shared the same MRCA may be more parsimoniously explained by a single ancestral introgression event between the branches that coalesce at this node, while those involving only a single species pair may have resulted from introgression between the extant species pair (Data S2). Another example of the former can be seen in clade 6 where the evidence suggests introgression occurred between the Hawaiian Scaptomyza and Drosophila (Figure S3) that are estimated to have diverged from each other more than 20 Mya. This ancient introgression may have occurred prior to the formation of Kauai island ~5 Mya which is now the oldest high island with extant species in these two groups55,56.
The matrix shows inferred introgression proportions as estimated from gene tree counts for the introgressed species pairs (STAR Methods), and then mapped to internal branches using the f-branch method20. The expanded tree at the top of each matrix shows both the terminal as well as ancestral branches. The tree on the left side of each matrix represents species relationships with mapped introgression events (red arrows) derived from the corresponding f-branch matrix (STAR Methods). The fractions next to each arrow represent the number of triplets that support a specific introgression event by both DCT and BLT divided by the total number of triplets that could have detected the introgression event. Dashed arrows represent introgression events with low support (triplet support ratio < 10%). See also Figures S3–S4, Tables S1–S2 and Data S2, S4–S5.
To summarize our DCT-BLT results and estimate both the number of introgression events and the proportion of the genome that introgressed during those events (![[Upsilon]](/corehtml/pmc/pmcents/x03D2.gif)
) we adapted the f-branch heuristic20 (implemented in Dsuite57; STAR Methods). Summed across all clades, our f-branch results suggest that at least 30 introgression events are required to explain our DCT-BLT results (Figure 3 and Figure S2). Clades 2, 4, 6, 7 and 9 showed the strongest evidence of introgression, in terms of both the total number of DCT-BLT significant triplets and estimates from Dsuite that support those events (Table 1). For example, in clade 2 Dsuite suggests an ancestral introgression event between the branch leading to D. obscura and D. bifasciata and the branch that leads to the clade containing D. pseudoobscura and D. affinis. Furthermore, this particular signal is characterized by a large fraction of introgressed genetic material ( = 0.237, Table 1) and by the large number of triplets that are significant according to both DCT and BLT (26 out of 28 total triplets that could detect this event are significant according to both tests). We stress that both our methods used to detect introgression (DCT and BLT) and our approaches for counting introgression events (f-branch) are conservative, and thus the true number of events could be substantially greater, as suggested by our analyses using QuIBL and HyDe. Regardless of the method used, careful examination of results in Data S2, Figure 3 and Figure S2 reveals that deep introgression events are clearly the best explanation for some of our patterns (e.g. the case from clade 7 involving D. robusta described above), although more recent events may have occurred as well (e.g. between D. pachea and D. acanthoptera; Data S2, clade 7).
Table 1.
Putative introgression events (![[long left and right arrow]](/corehtml/pmc/pmcents/x27F7.gif)
) are specified between different clades indicated by the pair of species in that clade with the oldest MRCA. The triplet ratio shows the number of significant and non-significant triplets according to DCT-BLT. The average γ was obtained from the f-branch results. The durations of the two introgressing lineages are represented by predicted lower and upper boundaries of credible intervals (95% CIs) estimated by MCMCTree using calibration scheme A.
| Clade | Introgression Event | Triplet ratio (significant/total) | Average γ | CI lower and upper bounds of lineage duration (Mya) |
|---|---|---|---|---|
| 2 | D. obscura…D. bifasciata D. pseudoobscura…D. affinis | 26/28 | 0.237 | 9.84–5.9, 11.53–6.18 |
| D. subobscura…D. guanche D. pseudoobscura…D. lowei | 2/40 | 0.016 | 9.84–1.93, 8.39–3.94 | |
| D. lowei D. azteca…D. affinis | 3/9 | 0.019 | 5.77–0, 8.39–2.88 | |
| 4 | D. ficusphila D. carrolli…D. elegans | 5/81 | 0.049 | 18.46–0, 16.72–8.55 |
| D. ficusphila D. erecta…D. eugracilis | 19/65 | 0.035 | 18.46–0, 14.77–10.3 | |
| D. erecta...D. orena D. mauritiana...D. melanogaster | 4/16 | 0.044 | 6.38–2.45, 7.71–2.92 | |
| 5 | D. leontia D. birchii...D. serrata | 1/55 | 0.076 | 3.07–0, 8.84–6.4* |
| 6 | S. pallida D. cyrtoloma...D. primaeva | 13/42 | 0.03 | 4.9–0, 24.22–4.52 |
| S. flava...S. montana D. cyrtoloma...D. prolacticillia | 25/56 | 0.032 | 16.14–1.87/6.56–3.05 | |
| D. primaeva D. cyrtoloma...D. silvestris | 1/40 | 0.02 | 6.56–0/3.96–2 | |
| D. heteroneura D. grimshawi...D. sproati | 1/36 | 0.021 | 1.1–0/3.22–2.17* | |
| D. primaeva D. prolacticillia | 1/12 | 0.012 | 6.56–0/2.16–0 | |
| 7 | D. robusta D. americana...D. montana | 49/60 | 0.113 | 29.36–0/26.88–4.28 |
| D. pseudotalamancana D. americana...D. montana | 52/55 | 0.103 | 23.56–0/26.88–4.28 | |
| D. novamexicana D. arizonae...D. hydei | 7/60 | 0.019 | 2.99–0/23.56–10.24 | |
| D. americana...D. novamexicana D. arizonae...D. seriema | 5/88 | 0.031 | 2.99–1.17/12.23–6.71* | |
| D. hydei D. americana...D. novamexicana | 2/28 | 0.034 | 13.82–0/2.99–1.17 | |
| D. hydei D. arizonae...D. seriema | 3/4 | 0.234 | 13.82–0/12.23–6.71 | |
| D. robusta D. pseudotalamancana | 6/16 | 0.076 | 29.36–0/23.56–0 | |
| D. pachea D. acanthoptera | 1/1 | 0.05 | 5.36–0/4.95–0 | |
| 8 | Z. camerounensis Z. lachaisei | 1/1 | 0.051 | 2.13–0/2.76–0 |
| Z. camerounensis Z. vittiger | 1/2 | 0.06 | 2.13–0/3–0 | |
| 9 | D. pruinosa D. arawakana...D. mush sãotomé | 47/110 | 0.138 | 22.41–0/27.26–18.21 |
| D. funebris...D. mush sãotomé D. albomicans…D. Pruinose | 206/288 | 0.031 | 22.7–14.74/27.26–16.86 | |
| D. subbadia D. guttifera…D. mush sãotomé | 1/18 | 0.106 | 3.19–0/19.51–11.12* | |
| D. innubila…D. mush sãotomé D. funebris | 1/2 | 0.135 | 15.94–7.48/19.51–0 | |
| D. immigrans D. neonasuta...D. sulfurigaster sulfurigaster | 1/24 | 0.01 | 1.44–0/3.5–1.7* | |
| D. immigrans (kari17) D. nasuta | 1/9 | 0.013 | 1.44–0/2.38–0 | |
| D. pallidipennis D. pulaua | 1/20 | 0.045 | 18.35–0/1.86–0 | |
| D. quadrilineata D. nasuta | 1/21 | 0.01 | 30.52–0/2.38–0 |
We note that some scenarios of ancestral population structure could potentially result in differences in the number and branch lengths of gene trees with either discordant topology (discussed in STAR Methods). We therefore applied a more stringent version of the DCT-BLT that compares the branch lengths of the discordant topology with those of the concordant topology; this test will not be sensitive to ancestral population structure but could potentially produce many false negatives (STAR Methods). The results of this test (Table S2) show that the majority (~2/3) of our strongly supported putative introgression events are inconsistent with the phenomenon of ancestral population structure-produced false positives. Given that this test is highly conservative, we interpret this result as evidence that the vast majority of our detected introgression events are true positives rather than artifacts of population structure.
To complement our f-branch analysis, we also used PhyloNet58,59 to identify the branches with the strongest signature of introgression in each of the nine monophyletic clades in our tree. Overall, PhyloNet’s results are largely consistent with our DCT-BCT analysis and therefore support our finding of multiple introgression events across the Drosophila phylogeny (see details in Figure S4).
Finally, we asked whether the proportion of the genome that introgressed between putatively introgressing taxa () varied with the timing of introgression events (Figure 4). Rather than timing introgression relative to when two hybridizing taxa shared a most recent common ancestor (which would require additional data, such as haplotype lengths of introgressed regions), we leveraged divergence time estimates across the drosophila phylogeny (Figure 1) and estimated when introgression events could have occured in time relative to the present (i.e., Mya). For this analysis, we focused on the 17 “best-supported” introgression events based on the criteria that more than 10% of the total triplets that could detect introgression between a given pair of taxa were significant according to both DCT and BLT (see solid red arrows in Figs. 3 and S3; Table 1). We estimated when these events occurred by taking the maximum, minimum, and midpoint times when the two branches that experienced introgression both coexisted in our dated phylogeny. We note that this approach results in imprecise time estimates, particularly for long branches in the phylogeny; however, it allowed us to test whether there was any obvious relationship between the proportion of the genome that introgressed () and when those introgression events took place in the past. In one instance, the two branches that putatively experienced introgression did not overlap in time in our phylogeny. This situation could be explained by “ghost” introgression with unsampled or extinct lineages. For the 17 remaining introgression events, there was not a significant relationship between the midpoint estimate of timing of introgression (Mya) and (Spearman’s rank correlation: 0.43; P = 0.085; Figure 4). Our analyses therefore support introgression across the evolutionary history of Drosophila, with introgressing species pairs exchanging a similar fraction of the genome (range of average estimates = 0.013 – 0.237) regardless of whether those events were ancient or more recent.
Each horizontal segment summarises one of the 17 introgression events highlighted in Figure 3 and is colored by clade. Segments span the times when the two putatively introgressing taxa both existed and are based on times inferred from the dating analysis summarised in Figure 1. Fraction of the genome that introgressed was estimated as the average f-branch statistic across all triplet comparisons that supported a given introgression event. Mya = million years ago.
Discussion
A time-calibrated tree of drosophilid evolution
Drosophila, as a genus, remains a premier model in genetics, ecology, and evolutionary biology. With over 1,600 species42, the genus has the potential to reveal why some groups are more speciose than others. Yet the phylogenetic relationships among the main groups in the genus have remained largely unresolved (reviewed in 42). Here we estimated a robust time-calibrated phylogeny for the whole genus using multilocus genomic data and calibrated it using a fossil record.
Our results confirm that the genus Drosophila is paraphyletic, with the genera Zaprionus, Scaptomyza, Leucophenga, and Hirtodrosophila each nested within the larger genus Drosophila. Consistent with the subdivisions previously proposed by refs. 60 and 44, clades 1–5 of our phylogeny contain species belonging to the subgenus Sophophora, and include species from the genus Lordiphosa (group A in Figure 1). Clades 6–9 of our phylogeny contain species belonging to the subgenus Drosophila (group B in Figure 1) and include species from the Hawaiian Drosophila and the subgenera Siphlodora, Phloridosa (synonymized with the subgenus Drosophila44, and genus Zaprionus. For more recent radiations within Drosophila, the topology we present is largely congruent with previous studies42,51 but two general observations are notable. First, our results confirm that Lordiphosa is closely related to the saltans and willistoni groups (clade 1) and part of the Sophophora subgenus (consistent with ref. 61). Second, we confirm that Zaprionus is related to the cardini/qunaria/immigrans group (consistent with refs. 42 and 60, but discordant with 43). Despite our well resolved phylogeny, comparisons with other studies emphasize the need to expand species sampling, especially given the potential to generate highly contiguous genomes at relatively low cost62.
Our results from divergence time analysis via MCMCTree suggest that the origin of Drosophila (including the subgenera Sophophora (group A) and Drosophila (group B)) occurred during the Eocene Epoch of the Paleogene, which is younger than estimates by 60 and 43, but older than estimates by 49. Different estimates of divergence times may be the result of different calibration information used, such as mutation rates, the time of formation of the Hawaiian Islands, and the fossil record. However, our comparison of various calibration schemes suggests that the choice of calibration information has a minor effect on MCMCTree’s age estimation (Figure S2). Additionally, credible intervals around our estimates tend to be notably narrower than in all of the aforementioned studies. In contrast to the previous studies, we used genome-scale multilocus data which would be expected to improve both the accuracy and precision of age estimates63,64.
On the other hand, we note that our analyses in BEAST2 using the FBD model yielded significantly older ages (Figure S2) especially for deeper nodes and with markedly wider credible intervals suggesting origination of Drosophila lineage in the Late Cretaceous. These calibration inconsistencies may arise as a result of the poor fossil record within Drosophila (only Scaptomyza dominicana from Dominican amber) and selection of the oldest fossils for deeper radiations, which together can lead to overestimation of nodal ages under the FBD model65. Moreover, the poor convergence behavior we observed would also be expected to produce larger credible intervals.
The extent of introgression in Drosophila
Access to genome-scale data has reinvigorated the study of hybridization and introgression14. We used genome-scale sequence data to provide the first systematic survey of introgression across the phylogeny of drosophilid flies. Our complementary—and conservative—approaches identified overlapping evidence for introgression within seven of the nine clades we analyzed (Figs. 3 and S3, Data S2). We conclude that at least 30 pairs of lineages have experienced introgression across Drosophila’s history (Table 1), though we note that other methods recover more introgression events (Data S2) and thus we cannot rule out the possibility that the true number is substantially higher. Moreover, we find that in many cases a substantial fraction of the genome is introgressed: our estimates indicate that numerous introgression events have altered gene tree topologies for >10% of the genome (Figs. 3 and S3, Table 1). Studies in contemporary Drosophila species suggest that selection may constrain the evolution of mixed ancestry, at least in naturally occurring9,36,66 and experimental admixed populations67,68. The results we have presented here used phylogenetic signals to show that introgression has nonetheless occurred and left a detectable signal within the genomes of many extant Drosophila.
In addition to providing an estimate of the extent of introgression, our results are informative about the timing of introgression among Drosophila lineages: the approaches we used to estimate the number of introgression events, and map them onto the phylogeny could potentially overestimate the timing of introgression if multiple independent more recent events are mistaken for one ancestral event. However, as described in the Results, both our PhyloNet analyses and a careful examination of our DCT-BLT results are most consistent with ancient introgression events in many cases. We also find evidence for very recent events, and although our analyses did not search for gene flow between sister taxa, previous studies of closely related species in Drosophila have revealed evidence of introgression9,10,29,31,32. Studies that have taken phylogenomic approaches to detect introgression in other taxa have also reported evidence for introgression between both “ancient” lineages (i.e., those that predate speciation events generating extant species) and extant species8,12,18,19,21. We conclude that introgression between Drosophila flies has similarly occurred throughout their evolutionary history.
Although the signal of introgression across our phylogeny provides evidence for widespread introgression in Drosophila, the evolutionary role of introgressed alleles remains to be tested. For example, the impact of hybridization and introgression on evolution can be diverse, from redistributing adaptive genetic variation23,69,70 to generating negative epistasis between alleles that have evolved in different genomic backgrounds (refs. 71–73; reviewed in refs. 15,16,74,75). The number of introgressed alleles that remain in a hybrid lineage depends on their selection coefficients76–78, their location in the genome (i.e., sex chromosomes vs. autosomes79–81), levels of divergence between the hybridizing species9,82,83, and recombination rates among loci6,84. Previous studies have, for example, shown that Drosophila hybrids often show maladaptive phenotypes36,85–89. Similarly, experimental hybrid swarms generated from two independent species pairs of Drosophila have shown that these populations can evolve to represent only one of their two parental species within as few as 10 generations, with the genome of one of their two parental species being rapidly purged from the populations67. These results show how hybrid Drosophila can be less fit than their parents, and further work is needed to determine the evolutionary effects, and the ecological context, of the introgression that we report here. However, our results suggest that not all introgressed material is deleterious in Drosophila, as we find that for some lineages a large fraction of the genome is introgressed (i.e. our estimates shown in Figs. 3 and S3 and Table 1). These results add to the growing body of literature that document a detectable phylogenetic signal of introgession left within the genomes of a wide range of species radiations that include Drosophila, other dipterans90, lepidopterans8,84,91, humans5,92,93, fungi1,2, and angiosperm plants11,12.
Caveats and future directions
We estimated the number of events required to explain the introgression patterns across the tree and in some cases those events were recovered as relatively ancient. However, our approaches for mapping gene flow events onto the phylogeny was somewhat parsimonious in that it favors older events over repeated and recent introgressions (see STAR Methods), and thus may bias the age of introgression towards ancient events and underestimate the true number of pairs of lineages that have exchanged genetic material. For example, introgression events we inferred at deeper nodes in our phylogeny are often supported by only a subset of comparisons between species pairs that spanned those nodes (e.g. see “ancient” introgression events in clades 2, 7 and 9; Figure 3). It is also possible that some patterns we observe reflect scenarios where introgressed segments have persisted along some lineages but been purged along others. This phenomenon could also cause older gene flow between sister lineages, which should generally be undetectable according to the BLT and DCT methods, to instead appear as introgression between non-sister lineages that our methods can detect. Future work could seek to more precisely reveal the number and timing of gene flow events across this phylogeny, including more recent introgression events and gene flow between extant and extinct/unsampled lineages, a pattern referred to as “ghost” introgression94,95.
Our analyses also do not identify the precise alleles that have crossed species boundaries or reveal the manner in which these alleles may have affected fitness in the recipient population74,75. Genome alignments, complete annotations, and/or population level sampling across the genus are required to determine whether certain genes or functional categories of genes are more likely to cross species boundaries than others. More complete taxonomic sampling, combined with methodological advances for inferring the number and timing of introgression events in large phylogenies, will increase our ability to identify the specific timing and consequences of introgression across Drosophila.
Conclusions
Speciation research has moved away from the debate of whether speciation can occur with gene flow to more quantitative tests of how much introgression occurs in nature, and how this introgression affects the fitness of individuals in the recipient population. Our well-resolved phylogeny and survey of introgression revealed that gene flow has been a relatively common feature across the evolutionary history of Drosophila. Yet, identifying the specific consequences of introgression on fitness and the evolution of species and entire radiations within Drosophila and other systems remains a major challenge. Future research could combine the power of phylogenomic inference with population-level sampling to detect segregating introgression between sister species to further our understanding of the amount, timing, and fitness consequences of admixture for diversification.
STAR METHODS
RESOURCE AVAILABILITY
Lead contact
More information regarding the resources used in this study should be directed to the lead contact Anton Suvorov (ude.cnu.dem@vorovus.ynotna)
Materials availability
This study did not generate any new unique reagents.
Data and code availability
The data and code produced during this study are publically available on GitHub (https://github.com/SchriderLab/drosophila_phylogeny) and FigShare (dx.doi.org/10.6084/m9.figshare.13264697). Whole genome sequencing data generated for this study are available on NCBI (BioProject PRJNA675888, BioProject PRJNA593822, and BioProject PRJNA611543).
METHOD DETAILS
Genome assemblies and public data
Genome sequences used by this work were obtained from concurrent projects and public databases. Genome sequencing and assembly for 84 genomes is described in62. These data are available for download at NCBI BioProject PRJNA675888. For the remaining genomes: sequencing and assembly of 8 Hawaiian Drosophila were provided by E. Armstrong and D. Price, described in Armstrong et al. (in prep) and available at NCBI BioProject PRJNA593822; sequences and/or assemblies of five nannoptera group species were provided by M. Lang and V. Courtier-Orgogozo and are available at NCBI BioProject PRJNA611543; 44 were downloaded as assembled sequences from NCBI GenBank; Z. sepsoides and D. neohypocausta were sequenced as paired-end 150bp reads on Illumina HiSeq 4000 at UNC and assembled using SPAdes v3.11.1 with default parameters96; and 15 were generated by assembling short read sequences downloaded from NCBI SRA. For sets of unassembled short reads, we used ABySS v2.2.397 with parameters “k=64” with paired-end reads (typically 100–150bp) to assemble the reads. Finally, outgroup genome sequences (A. gambiae, M. domestica, L. trifolii, C. hians, and E. gracilis) were obtained from NCBI GenBank. See Data S3 for a full list of samples, strain information, accessions, and associated publications.
Orthology Inference
We identified single-copy orthologous genes in each genome using BUSCO (Benchmarking Universal Single-Copy Orthologs; v3.1.098). BUSCO was run with orthologs from the Diptera set in OrthoDB v.9 (odb9) using default parameters. For each species, all BUSCOs found in a single copy were used for phylogenetic analysis.
Assignment of BUSCO genes to Muller elements for obscura group species
Each of the BUSCO genes identified as single-copy in each of the group 12 (obscura group: D. affinis, D. athabasca, D. azteca, D. bifasciata, D. guanche, D. lowei, D. miranda, D. obscura, D. persimilis, D. pseudoobscura, D. subobscura) genome assemblies was assigned to one of the six Muller elements (elements A-F). For D. athabasca, D. bifasciata, D. lowei, D. miranda, D. pseudoobscura, and D. subobscura, contig/scaffold associations with chromosomes and/or Muller elements were simply obtained from NCBI GenBank assembly report tables. For the remaining genomes (D. affinis, D. azteca, D. guanche, D. obscura, D. persimilis), we used whole-genome alignments to infer the Muller element associated with each contig or scaffold. Using the Progressive Cactus99 software, each remaining genome was aligned to a closely related reference genome (D. affinis - D. athabasca; D. azteca - D. athabasca; D. guanche - D. subobscura; D. obscura - D. bifasciata; D. persimilis - D. miranda) with a similar karyotype54,100. Using the reference genomes as backbones, each remaining genome was then scaffolded, with Ragout101. The scaffolds allowed us to annotate each contig in the remaining genomes with Muller element information from the reference genomes (see Data S4). BUSCO genes on unplaced contigs were ignored.
Phylogenetic reconstruction
Every DNA BUSCO locus was aligned with MAFFT v7.427102 using the L-INS-i method. We removed sites that had fewer than three non-gap characters from the resulting multiple sequence alignments (MSAs). These trimmed MSAs were concatenated to form a supermatrix. To assess the quality of the assembled supermatrices we computed pairwise completeness scores in AliStat103 (Figure S5). We inferred a maximum likelihood (ML) phylogenetic tree from the supermatrix (a.k.a. concatenated alignment) using IQ-TREE v1.6.5104, and treated the supermatrix as a single partition. IQ-TREE was run under GTR+I+G substitution model, as inference under any other substitution model will not necessarily lead to better accuracy of tree topology estimation105. To estimate the support for each node in this tree, we used three different reliability measures. We did 1,000 ultrafast bootstrap (UFBoot) replicates106 and additionally performed an approximate likelihood ratio test with the nonparametric Shimodaira–Hasegawa correction (SH-aLRT) and a Bayesian-like transformation of aLRT107. We used the ML gene trees obtained by IQ-TREE with a GTR+I+G substitution model for tree inference in ASTRAL108. For the estimated ASTRAL tree we calculated the support of each node using local posterior probabilities (LPP)108. Also, we created a gene tree set by removing taxa with outlier branch lengths that were potentially produced by misaligned regions and/or incorrect orthology inference in TreeShrink45 under default parameters. This analysis resulted in a small fraction of branches removed from our gene tree set (<5.5%)
We did two additional analyses to verify the robustness of our topology inference. First, we inferred an ML tree using WAG+I+G substitution model from the protein supermatrix obtained from concatenation of protein BUSCO MSAs. MSAs based on amino acid sequences have been shown to have superior accuracy to DNA MSAs for distantly related species109. Second, to verify that long branch attraction did not distort our tree topology, we inferred an ML tree under a GTR+I+G substitution model using a different set of outgroup species from the DNA supermatrix. Specifically, instead of distantly related Anopheles gambiae, we used Musca domestica, Liriomyza trifolii, Curricula hians and Ephydra gracilis together as our outgroup species.
Phylogenetic Support Analysis via Quartet Sampling
We used quartet sampling (QS) as an additional approach to estimate phylogenetic support11. Briefly, QS provides three scores for internal nodes: (i) quartet concordance (QC), which gives an estimate of how sampled quartet topologies agree with the putative species tree; (ii) quartet differential (QD) which estimates frequency skewness of the discordant quartet topologies, and can be indicative of introgression if a skewed frequency observed, and (iii) quartet informativeness (QI) which quantifies how informative sampled quartets are by comparing likelihood scores of alternative quartet topologies. Finally, QS provides a score for terminal nodes, quartet fidelity (QF), which measures a taxon “rogueness”. We did QS analysis using the DNA BUSCO supermatrix described above, specifying an IQ-TREE engine for quartet likelihood calculations with 100 replicates (i.e. number of quartet draws per focal branch).
Fossil Dating
MCMCTREE:
We implemented the Bayesian algorithm of MCMCTree v4.9h46 with approximate likelihood computation to estimate divergence times within Drosophila using several calibration schemes (Data S1). First, we estimated branch lengths by ML and then the gradient and Hessian matrices around these ML estimates in MCMCTree using the DNA supermatrix and species tree topology estimated by IQ-TREE. Because large amounts of sequence data are not essential for accurate fossil calibration110, we performed dating analysis using a random sample of 1,000 MSA loci (out of 2,791) for the sake of computational efficiency. Thus, for this analysis the supermatrix was generated by concatenating 1,000 randomly selected gene-specific MSAs. Using fewer loci (10 and 100) for fossil calibration did not drastically affect nodal age estimation (Figure S1). We removed sites that had less than 80 non-gap characters from all these supermatrices. Second, we used the gradient and Hessian matrix, which constructs an approximate likelihood function by Taylor expansion111, to perform fossil calibration in MCMC framework. For this step we specified a GTR+G substitution model with four gamma categories; birth, death and sampling parameters of 1, 1 and 0.1, respectively. To model rate variation we used an uncorrelated relaxed clock. To ensure convergence, the analysis was run ten times independently for 8 × 106 generations (first 106 generations were discarded as burn-in), logging every 1,000 generations. We used the R package MCMCtreeR112 to visualize the calibrated tree.
BEAST 2:
Additionally we performed fossil calibration using the Fossilized Birth-Death (FBD) process47 as implemented in the Bayesian framework of BEAST 2.6.348. For scalability purposes, we randomly selected 1,000 loci and then partitioned them into 10 supermatrices each consistent of 100 different MSAs. Each of these 10 datasets was treated as a single partition in the downstream analyses. Additionally, we removed sites that had less than 128 non-gap characters from all these supermatrices. To perform fossil calibration, we used a GTR+G model with four gmamma categories, and an optimized relaxed clock113 was used to model rate variation. For the FBD prior we specified an initial origin value of 230 Mya (which corresponds to the age of oldest known dipteran fossil Grauvogelia), and the tree likelihood was conditioned on the proportion of species sampled at present (p= 0.1). The remaining priors were set to their defaults. In order to directly compare divergence time estimation between BEAST 2 and MCMCTree, we used the same fixed IQ-TREE species tree topology with several exceptions. First, we did not fix the phylogenetic positions of contemporary Scaptomyza species and fossil taxon Scaptomyza dominicana within its monophyletic group. Second, we did not constrain relationships of outgroup species L. varia, C. costata, S. lebanonensis including fossil taxon Electrophortica succini. Two additional fossils, Oligophryne and Phytomyzites, were specified for Drosophilidae stem. Furthermore, to accomodate uncertainty of fossil dates we incorporated age ranges for several fossils (Data S1). For each of the 10 datasets we ran 2 independent MCMC chains for 6 × 108 generations with sampling frequency of 10,000 for each model parameter. Additionally, we performed sampling from the prior distribution only. Convergence was assessed using ESS in Tracer114. Divergence times were generated by taking means of posterior nodal ages discarding 25% of the sampled trees as burn-in in TreeAnnotator for each dataset. To drastically improve computational efficiency of likelihood calculations in all BEAST 2 analyses we used the program in conjunction with BEAGLE library115 that enables GPU utilization.
Inferring Introgression Across the Tree
Gene tree-based methods:
In order to detect patterns of introgression we used three different methods that rely on the topologies of gene trees, and the distributions of their corresponding branch lengths, for triplets of species. If the true species tree is ((A, B), C), these tests are able to detect cases of introgression between A and C, or between B and C. These include two of the methods that we devised for this study, and which use complementary pieces of information—the counts of loci supporting either discordant topology, and the branch-length distributions of gene trees supporting these topologies, respectively—to test an introgression-free null model.
The first method we developed was the discordant-count test (DCT), which compares the number of genes supporting each of the two possible discordant gene trees: ((A, C), B) or (A, (B, C)), similar in principle to the delta statistic from116. Genes may support the two discordant topologies (denoted T1 and T2) in the presence of ILS and/or in the presence of introgression. In the absence of ancestral population structure, gene genealogies from loci experiencing ILS will show either topology with equal probability; ILS alone is not expected to bias the count towards one of the topologies. In the presence of introgression, one of the two topologies will be more frequent than the other because the pair of species experiencing gene flow will be sister lineages at all introgressed loci (illustrated in Figure 2). For example, if there is introgression between A and C, there will be an excess of gene trees with the ((A, C), B) topology. The DCT identifies pairs of species that may have experienced introgression by performing a χ2 goodness-of-fit test on the gene tree count values for a species triplet to determine whether their proportions significantly deviate from 0.5, the expected proportion for each gene genealogy under ILS. We used this test on all triplets extracted from BUSCO gene trees within each clade, and the resulting P-values were then corrected for multiple testing using the Benjamini-Hochberg procedure with a false discovery rate (FDR) cutoff of 0.05. We note that these tests are not independent since different triplets may contain overlapping taxa. Thus, while our correction results in more conservative tests57, the inferred FDRs may be somewhat inaccurate.
Second, we devised a branch-length test (BLT) to identify cases of introgression (illustrated in Figure 2). This test examines branch lengths to estimate the age of the most recent coalescence event (measured in substitutions per site). Introgression should result in more recent coalescences than expected under the concordant topology with complete lineage sorting, while ILS shows older coalescence events90. Importantly, ILS alone is not expected to result in different coalescence times between the two discordant topologies, and this forms the null hypothesis for the BLT. For a given triplet, for each gene tree we calculated the distance d (a proxy for the divergence time between sister taxa) by averaging the external branch lengths leading to the two sister taxa under that gene tree topology. We calculated d for each gene tree and denote values of d from the first discordant topology dT1 and those from the second discordant topology dT2. We then compared the distributions of dT1 and dT2 using a Mann-Whitney U test. Under ILS alone the expectation is that dT1 = dT2, while in the presence of introgression dT1 < dT2 (suggesting introgression consistent with discordant topology T1) or dT1 > dT2 (suggesting introgression with consistent with topology discordant T2). The BLT is conceptually similar to the D3 test117, which transforms the values of dT1 and dT2 in a manner similar to the D statistic for detecting introgression92. As with the DCT, we performed the BLT on all triplets within a clade and used a Benjamini-Hochberg correction with a false discovery rate cutoff (FDR) of 0.05. We note that both the DCT and BLT will be conservative in cases where, for a triplet ((A,B), C), there is introgression between A and C as well as B and C, with the extreme case of equal rates of introgression for both species pairs resulting in a complete loss of power.
Finally, we used QuIBL8, an analysis of branch-length distribution across gene trees to infer putative introgression patterns. Briefly, under coalescent theory internal branches of rooted gene trees for a set of 3 taxa (triplet) can be viewed as a mixture of two distributions: one that generates branch lengths under ILS, and the other under introgression/speciation. Thus, the estimated mixing proportions (π1 for ILS and π2 for introgression/speciation; π1 + π2 = 1) of those distribution components show which fraction of the gene trees were generated through ILS or non-ILS processes. For a given triplet, QuIBL computes the proportion of gene trees that support the three alternative topologies. Then for every alternative topology QuIBL estimates mixing proportions along with other relevant parameters via Expectation-Maximization and computes Bayesian Information Criterion (BIC) scores for ILS-only and introgression models. For concordant topologies elevated values of π2 are expected whereas for discordant ones π2 values significantly greater than zero are indicative of introgression. To identify significant cases of introgression here we used a cutoff of ΔBIC < −30 as in8. We ran QuIBL on every triplet individually under default parameters with the number of steps (the numsteps parameter) set to 50 and using Anopheles gambiae for triplet rooting; the branch length between A. gambiae and the triplet is not used for any of QuIBL’s calculations.
We note that the DCT and BLT methods are potentially impacted by ancestral population structure: if the lineages leading to B and C were in subpopulations that were more likely to interbreed in the ancestral population, then the ((B, C), A) topology might be expected to be more prevalent than ((A, C), B), along with a shorter time back to the first coalescence. However, it is unclear how much of a concern ancestral population structure should be for this analysis, as it seems less likely that it would be a pair of lineages that diverged first (i.e. A and C or B and C) that interbred more frequently in the ancestral population instead of the two lineages that went on to be sister taxa (i.e. A and B). Nonetheless, plausible scenarios of ancestral structure supporting one discordant topology over the other can be devised (e.g. ref 118). We therefore conducted a more stringent version of our DCT-BLT combined test that requires the average distance between the two introgressing taxa (when examining gene trees with the discordant topology consistent with introgression) to be less than that between the two sister species (when examining gene trees with the concordant topology). Such a pattern is consistent with introgression between non-sister species, which must occur more recently than the species split and therefore causing more recent coalescence events, but not with ancestral structure which will still result in older coalescence times for discordant trees than the concordant trees (because structure in the ancestral population is only a factor in the case of ILS). Note that this test is expected to be especially conservative because ILS, which for many triplets accounts for a sizable fraction of our discordant gene trees, will push the coalescent times for all discordant topologies back further in time.
We also examined the effect of evolutionary rate heterogeneity measured in branch-specific dN/dS values on introgression detection. To that end, we generated codon alignments for each BUSCO locus using TranslatorX119 and then calculated dN/dS ratios for each gene tree in PAML46 within each clade using a free-ratios branch model that assumes independent dN/dS for each gene tree branch. Then, we evaluated the distribution of dN/dS ratios across all gene trees to determine the 95th percentile value of dN/dS. Thus, we repeated our DCT/BLT analyses for each triplet after excluding every gene tree that had at least one branch with dN/dS > 0.53. Note, branches with dN/dS values where dS< 0.001 or >5 were deemed unreliable and thus were excluded from calculation of a critical value or from downstream filtering. Additionally, we performed random filtering of gene trees to see if this procedure would have a similar impact on downstream introgression-detection as did our dN/dS filter. First, we estimated the distribution of proportions of gene trees retained for each triplet after applying the dN/dS filter. Then, for a given triplet, we randomly drew a number of genes to remove from the aforementioned distribution, and then applied our DCT-BLT method to this triplet after removing the selected number of genes. This process was repeated for each triplet tested in our main analysis to generate a randomly filtered set of DCT-BLT results for each of our 9 clades. We then repeated this entire process 1000 times and noted the average fraction of DCT-BLT results remaining significant after randomly filtering genes.
Our DCT-BLT test assumes that there is no recombination within loci and complete inter-locus independence—these assumptions are commonly made by introgression inference methods10,120,121. We note that intra-locus recombination may interfere with the signatures of introgression by reducing discordant topology counts (because even loci experiencing introgression will have non-introgressed segments), and similarly diluting branch-length signatures of introgression, thereby reducing the sensitivity of our DCT-BLT approach. Nevertheless, site-pattern-based approaches (e.g. HyDe, see below) are not affected by intra-locus recombination as they evaluate each site in an MSA independently.
Site-pattern -based detection of introgression:
Signatures of introgression can be identified by investigating fractions of certain site patterns within MSAs of species quartets. One of the most widely used methods is based on the counts of ABBA-BABA site patterns (aka., Patterson’s D statistic122). Here we used the hybridization model implemented in HyDe52 that implements an alternative invariant-based statistic to test introgression and estimate the fraction of the introgressed genome (). We ran HyDe analysis on each of the 9 clades using the entire supermatrix and in each case selected the quartet’s outgroup from a sister clade. Additionally, to examine effects of outgroup choice, we ran HyDe analyses with a more distantly related outgroup, Anopheles gambiae for all clades. The resulting P-values for each quartet were corrected for multiple testing using the Bonferroni method. To investigate an individual contribution of each BUSCO locus to introgression, we additionally ran HyDe using BUSCO MSAs with Anopheles gambiae outgroup. We note, however, in this case HyDe’s power to detect introgression will be reduced, especially for short MSAs with <10,000 sites52. A complete summary for each BUSCO locus including introgression results from locus-specific HyDe and BLT/DCT analyses is included in Data S5.
Placing introgression events on the phylogeny:
All the aforementioned methods can infer multiple correlated signatures of introgression especially when triplets/quartets share the same taxa. Thus it can be difficult to interpret these interdependent results. To alleviate this problem,20 devised a simple heuristic metric called f-branch to disentangle and map introgression events detected in multiple correlated species pairs onto the tree. In the original formulation, f-branch examines multiple f4 statistics measured for each species pair and that quantify, the proportion of introgressed material for that pair. However, the calculation of the f4 statistic requires allele frequency measures within each sampled species. Thus, to calculate f-branch statistic, instead of f4 we used the introgression proportion derived from DCT/BLT as follows: estimates and to visualize the results within each clade we used the Dsuite python package57.
Dsuite outputs a matrix of estimates that have been partially collapsed: on one axis of this matrix signals of introgression can appear on ancestral branches, but on the other axis only extant branches are shown. Thus, we manually further collapsed these signatures by parsimoniously assuming that if some lineage A showed evidence of introgression with multiple descendants of some other lineage B that is not ancestral to A, then we considered this to be caused by a single introgression event between A and B. Note that we did not require all descendants of lineage B to share this signature of introgression, and thus this approach could potentially undercount the number of introgression events and overestimate their ages.
Phylogenetic networks:
Introgression generates instances of reticulate evolution such that purely bifurcating trees cannot adequately represent evolutionary history; phylogenetic networks have been shown to provide a better fit to describe these patterns123,124. We used PhyloNet58,59 to calculate likelihood scores for networks generated by placing a single reticulation event (node) in an exhaustive manner, i.e. connecting all possible branch pairs within a clade and determining which of the resulting phylogenetic networks produced the best likelihood score. We note that networks with more reticulation events would most likely exhibit a better fit to observed patterns of introgression but the biological interpretation of complex networks with multiple reticulations is more challenging; thus, we limited the analysis to a single reticulation event even though this will produce false negatives in clades with multiple gene flow events. Because full likelihood calculations with PhyloNet can be prohibitively slow for large networks, for each of clades 1 through 9 we selected a subsample of 10 species in a manner that preserves the overall species tree topology. No subsampling was performed for clade 3 which has fewer than 10 species. Using these subsampled clade topologies, we formed all possible network topologies having a single reticulation node (with the exception of networks having reticulation nodes connecting sister taxa). Because PhyloNet takes gene trees as input, for each clade we subsampled each gene tree to include only the subset of 10 species selected for the PhyloNet analysis (or all species in the case of clade 3); any gene trees missing at least one of these species were omitted from the analysis. Finally, we used the GalGTProb program125 of the PhyloNet suite to obtain a likelihood score for each network topology for each clade. We report networks with the highest likelihood scores.
QUANTIFICATION AND STATISTICAL ANALYSIS
Nodal reliability measures including ultrafast bootstrap (UFBoot), an approximate likelihood ratio test with the nonparametric Shimodaira–Hasegawa correction (SH-aLRT) and a Bayesian-like transformation of aLRT were used to assess phylogenetic nodal support. The DCT approach conducts a χ2 goodness-of-fit test for an imbalance in the number of discordant gene trees, and the BLT approach uses the Mann-Whitney U test to assess the significance of differences in coalescent times between the discordant gene trees (see “STAR Methods” for a detailed description of these tests). Multiple testing correction techniques of a Benjamini-Hochberg with a false discovery rate (FDR) and Bonferroni were used to adjust P-values derived from DCT/BLT and HyDe, respectively, and a significance threshold of 0.05 was applied to these adjusted P-values. The difference of Bayesian Information Criterion scores (ΔBIC) was used to assess significance of QuIBL results at a cutoff of −30. The correlation between timing of introgression events and was tested using Spearman’s rank correlation.
KEY RESOURCES TABLE
Supplementary Material
2
Data S2. Results of introgression analyses (DCT, BLT, DCT-BLT, HyDe and QuIBL) for tested species pairs, related to Figure 3.Columns show tested species pairs, number of triplets that exhibit significant introgression signal for species pair, total number of tested triplets and fraction of significant vs. total tested triplets.
3
Data S3. Taxon sampling and genome assembly information, related to Figure 1.
4
Data S5. Correspondence between BUSCO loci and Muller elements, related to Figure 3.
7
Acknowledgements
We thank M. Hahn, M. Turelli, A. Yassin, D. Obbard and M. Matschiner for helpful feedback on a previous draft, and M. Hibbins for sharing simulated gene trees. AS and DRS were supported by the NIH under award nos. R00HG008696 and R35GM138286. BYK was supported by the NIH under award no. F32GM135998. JW was supported by the NIH under award no. K01DK119582. DP, AAC were supported by NSF Dimensions of Biodiversity award 1737752. DRM and ERRD were supported by NIH award R01GM121750. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Footnotes
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Declaration of interests
The authors declare that there is no conflict of interest
REFERENCES
Citations & impact
Impact metrics
Citations of article over time
Alternative metrics
Discover the attention surrounding your research
https://www.altmetric.com/details/117014853
Smart citations by scite.ai
Explore citation contexts and check if this article has been
supported or disputed.
https://scite.ai/reports/10.1016/j.cub.2021.10.052
Article citations
Extreme genome scrambling in marine planktonic Oikopleura dioica cryptic species.
Genome Res, 34(3):426-440, 25 Apr 2024
Cited by: 0 articles | PMID: 38621828 | PMCID: PMC11067885
Hybrid speciation driven by multilocus introgression of ecological traits.
Nature, 628(8009):811-817, 17 Apr 2024
Cited by: 1 article | PMID: 38632397 | PMCID: PMC11041799
Horizontal transfer and the widespread presence of Galileo transposons in Drosophilidae (Insecta: Diptera).
Genet Mol Biol, 46(3 suppl 1):e20230143, 29 Mar 2024
Cited by: 0 articles | PMID: 38569056 | PMCID: PMC10990002
A comprehensive evolutionary scenario for the origin and neofunctionalization of the Drosophila speciation gene Odysseus (OdsH).
G3 (Bethesda), 14(3):jkad299, 01 Mar 2024
Cited by: 0 articles | PMID: 38156703 | PMCID: PMC10917504
A new hybrid incompatibility locus between Drosophila melanogaster and Drosophila sechellia.
Genetics, 226(3):iyae001, 01 Mar 2024
Cited by: 1 article | PMID: 38184848
Go to all (62) article citations
Data
Data behind the article
This data has been text mined from the article, or deposited into data resources.
Data Citations
- (1 citation) DOI - 10.6084/m9.figshare.13264697
BioProject (3)
- (2 citations) BioProject - PRJNA675888
- (2 citations) BioProject - PRJNA593822
- (2 citations) BioProject - PRJNA611543
Similar Articles
To arrive at the top five similar articles we use a word-weighted algorithm to compare words from the Title and Abstract of each citation.
Primate phylogenomics uncovers multiple rapid radiations and ancient interspecific introgression.
PLoS Biol, 18(12):e3000954, 03 Dec 2020
Cited by: 38 articles | PMID: 33270638 | PMCID: PMC7738166
Phylogenomic approaches untangle early divergences and complex diversifications of the olive plant family.
BMC Biol, 20(1):92, 25 Apr 2022
Cited by: 22 articles | PMID: 35468824 | PMCID: PMC9040247
Reticulate evolutionary history and extensive introgression in mosquito species revealed by phylogenetic network analysis.
Mol Ecol, 25(11):2361-2372, 10 Mar 2016
Cited by: 54 articles | PMID: 26808290 | PMCID: PMC5146328
Incongruence in the phylogenomics era.
Nat Rev Genet, 24(12):834-850, 27 Jun 2023
Cited by: 13 articles | PMID: 37369847
Review
Phylogenomic approaches to detecting and characterizing introgression.
Genetics, 220(2):iyab173, 01 Feb 2022
Cited by: 42 articles | PMID: 34788444 | PMCID: PMC9208645
Review
Free full text in Europe PMC
Funding
Funders who supported this work.
NHGRI NIH HHS (1)
Grant ID: R00 HG008696
NIDDK NIH HHS (1)
Grant ID: K01 DK119582
NIGMS NIH HHS (5)
Grant ID: R35 GM138286
Grant ID: F32 GM135998
Grant ID: R35 GM118165
Grant ID: R01 GM121750
Grant ID: T32 GM128596
NSF
National Institutes of Health (5)
Grant ID: K01DK119582
Grant ID: F32GM135998
Grant ID: R35GM138286
Grant ID: R00HG008696
Grant ID: R01GM121750
National Science Foundation (1)
Grant ID: 1737752