Phylotranscriptomics Reveals Discordance in the Phylogeny of Hawaiian Drosophila and Scaptomyza (Diptera: Drosophilidae).

Island radiations present natural laboratories for studying the evolutionary process. The Hawaiian Drosophilidae are one such radiation, with nearly 600 described species and substantial morphological and ecological diversification. These species are largely divided into a few major clades, but the relationship between clades remains uncertain. Here, we present new assembled transcriptomes from 12 species across these clades, and use these transcriptomes to resolve the base of the evolutionary radiation. We recover a new hypothesis for the relationship between clades, and demonstrate its support over previously published hypotheses. We then use the evolutionary radiation to explore dynamics of concordance in phylogenetic support, by analyzing the gene and site concordance factors for every possible topological combination of major groups. We show that high bootstrap values mask low evolutionary concordance, and we demonstrate that the most likely topology is distinct from the topology with the highest support across gene trees and from the topology with highest support across sites. We then combine all previously published genetic data for the group to estimate a time-calibrated tree for over 300 species of drosophilids. Finally, we digitize dozens of published Hawaiian Drosophilidae descriptions, and use this to pinpoint probable evolutionary shifts in reproductive ecology as well as body, wing, and egg size. We show that by examining the entire landscape of tree and trait space, we can gain a more complete understanding of how evolutionary dynamics play out across an island radiation.

Introduction

In the era of genome-scale data, we have an opportunity to unpack the biological meaning of phylogenetic support. In phylogenetic analyses that seek to discover the relationships between organisms, support is often defined as the proportion of information that favors a particular branch in an evolutionary tree (Simon 2020). Methods have been developed that emphasize extracting the tree with the greatest amount of support from out of an otherwise rugged landscape of treespace (Swofford et al. 1996; Hedtke et al. 2006). However, a growing number of studies have emphasized the biological relevance of that landscape to our understanding of the evolutionary process (Kellogg et al. 1996; Baum 2007; Smith et al. 2015). For example, many new studies have contributed evidence that, even with trees with high measures of conventional support, we can expect large amounts of discordance among sites and genes, especially when examining speciation events with short internodes or with a likelihood of introgression (Weisrock et al. 2012; Pease et al. 2016). Here, we use the island radiation of Hawaiian drosophilid flies to study the landscape of treespace, and show that the relationships between the major groups of these flies are best understood by using methods that embrace evolutionary discordance.

The Hawaiian Drosophila have a long history as a model clade for the implementation of phylogenetic methods (O’Grady and DeSalle 2018). More than 20 years ago, Baker and DeSalle (1997) used the Hawaiian radiation of Drosophila to perform one of the first analyses to demonstrate incongruence between an overall species tree and underlying gene trees. Their study focused on the resolution between major clades of Hawaiian Drosophila and built on the landmark work done by Carson in the 1970s inferring the phylogeny of a subgroup of Hawaiian Drosophila, the picture-wing flies, based on the banding pattern of polytene chromosomes (Carson and Kaneshiro 1976), among other early phylogenetic studies (Beverley and Wilson 1985; Thomas and Hunt 1991). During the past 20 years, the relationships between major groups have been revisited several times (Kambysellis et al. 1995; Bonacum 2001). Most recently, O’Grady et al. (2011) used mitochondrial genes and expanded taxon sampling, and Magnacca and Price (2015) used an expanded nuclear gene set. The study presented here builds on this foundational work, presenting the first phylogenetic analysis of genome-scale data for the group.

The Hawaiian Drosophilidae consist of 566 described species (O’Grady et al. 2010; Magnacca and Price 2012), with hundreds more estimated to be awaiting description (O’Grady et al. 2010). These species have been divided into the following major clades (O’Grady et al. 2010): 1) the picture-wing, nudidrosophila, ateledrosophila (PNA) clade, which has served as a model clade for the study of sexual selection (Kaneshiro and Boake 1987) and speciation (Kang et al. 2016); 2) the antopocerus, modified-tarsus, ciliated-tarsus (AMC) clade, first proposed by Heed (1968) and (O’Grady et al. 2010) and confirmed by subsequent phylogenetic studies (O’Grady et al. 2011; Lapoint et al. 2014); 3) the modified-mouthparts (MM) clade; and 4) the haleakalae clade, an enigmatic group in need of further study (Hardy et al. 2001). Several other smaller clades have been suggested as falling outside of these major groups, including the rustica group of three species (O’Grady et al. 2001), and the monotypic lineages of D. primaeva and D. adventitia. The position of D. primaeva has been somewhat uncertain, but several studies have suggested it is the sister taxon to picture-wing flies (Bonacum 2001), including the work on polytene chromosomes by Carson and Stalker (1969). The species D. adventitia was originally suggested to be part of the MM clade (Hardy 1965), but recent studies placed it as the sister taxon to D. primaeva (Bonacum 2001) or possibly other major clades. Additionally, the Hawaiian Drosophila are the sister clade of the genus Scaptomyza, which is nested within the broader paraphyletic genus Drosophila and is hypothesized to have colonized the island independently (Throckmorton 1966; Lapoint et al. 2013), possibly more than once (Katoh et al. 2017). Throughout this manuscript, we use Hawaiian Drosophila to refer to non-Scaptomyza Hawaiian species, and Hawaiian Drosophilidae to refer to the clade of Hawaiian Drosophila+Scaptomyza.

Many phylogenetic studies have been performed which have confirmed the monophyly of each of these clades and provided resolution for internal relationships (PNA, Bonacum et al. 2005; Magnacca and Price 2015; AMC, Lapoint et al. 2011, 2014; haleakalae, O’Grady and Zilversmit 2004; and Scaptomyza, Lapoint et al. 2013; Katoh et al. 2017). Previous phylogenetic studies, however, have not resulted in a consensus relationship between the major clades within Hawaiian Drosophila (supplementary fig. S1, Supplementary Material online) (Magnacca and Price 2015). Magnacca and Price (2015) showed that different phylogenetic methods of analysis (e.g., using software based on Bayesian statistics rather than maximum likelihood for inference) produced highly incongruent topologies (fig. 1) (Magnacca and Price 2015). In that study, the most likely topology had D. primaeva as the sister taxon to all other Hawaiian Drosophila, and included a clade uniting MM+AMC+haleakalae, with the haleakalae clade showing greater affinity to AMC species relative to MM species (fig. 1). This topology was consistent with the tree suggested by O’Grady et al. (2011) analyzing mitochondrial data and using maximum likelihood and Bayesian analyses. However, the analyses of Magnacca and Price (2015) using the Bayesian software package BEAST showed an alternative relationship, with haleakalae flies as the sister clade to all other Hawaiian Drosophila, a clade uniting the MM+PNA+D. primaeva, and closer affinity between D. primaeva and PNA species than between D. primaeva and MM species (fig. 1). This latter arrangement is largely consistent with relationships proposed by Throckmorton (1966) and reiterated in several subsequent studies (supplementary fig. S1, Supplementary Material online) (Kambysellis et al. 1995; Baker and DeSalle 1997; Bonacum 2001).

Fig. 1.

Phylotranscriptomic analysis indicates relationships between major clades distinct from those previously hypothesized. Photos show six of the 12 species with de novo transcriptomes presented in this study, listing their parent clade and the Hawaiian island on which they are found. (A) Results novel to this study, showing best supported tree across maximum likelihood and Bayesian analyses. Node labels indicate ultrafast bootstrap values, gene tree concordance factors (gCF), and site concordance factors (sCF), see concordance factor analysis below. Drosophila adventitia was not present in phylotranscriptomic analyses; see figure 3 for information on its placement. (B, C) Previously hypothesized relationships between the picture wing-nudidrosophila-ateledrosophila (PNA), modified-mouthparts (MM), antopocerus-modified tarsus-ciliated tarsus (AMC), haleakalae, and Scaptomyza clades, as well as two monotypic clades, D. primaeva and D. adventitia. Topology B was recovered in O’Grady et al. (2011) and Magnacca and Price (2015). Topology C was recovered using the Bayesian software BEAST in Magnacca and Price (2015), showing incongruent relationships between clades at the base of the radiation of Hawaiian Drosophila.

Fig. 3.

Time-calibrating the phylogeny of 316 Drosophilidae species. This phylogeny was inferred using IQtree to analyze all publicly available genetic data for Hawaiian Drosophila and Scaptomyza. It was then calibrated to time using the software BEAST, with four calibration points at nodes that show a biogeographic progression rule (Magnacca and Price 2015). Similar results were obtained using a calibration scheme (Russo et al. 2013) that takes into account a single fossil taxon for the group (table 1). The 95% highest posterior density intervals for each node are shown as gray bars, indicating the credible interval for the age of that group. The age at which four Hawaiian islands are estimated to have become habitable is shown in green. Colored labels indicate the clade to which taxa belong, and colors correspond to figure 1; taxa without a colored label are species with genetic data that are as of yet undescribed. See supplementary figure S13, Supplementary Material online, for bootstrap support. Calibration using only island biogeography is known to be imprecise (Kodandaramaiah 2011); therefore, the divergence times shown here are considered tentative.

Resolving these relationships is critical for our understanding of the morphological and ecological evolution of these flies (Kambysellis et al. 1995; Bonacum 2001; O’Grady et al. 2011). Hawaiian Drosophila demonstrate a large diversity in body size (Stevenson et al. 1995), wing size (Edwards et al. 2007), and egg size (Montague et al. 1981); in the number and position of structural features such as wing spots (Edwards et al. 2007); in the number of egg-producing units in the ovary (ovarioles) (Kambysellis and Heed 1971; Sarikaya et al. 2019); and in the type of substrate used for oviposition and larval feeding (Kambysellis et al. 1995; Magnacca et al. 2008). Some clades demonstrate unique suites of morphological and behavioral traits, whose evolutionary history is unclear because of uncertainties in the phylogeny. For example, the haleakalae flies exclusively use fungal oviposition substrates and are considered to have less complex mating behaviors than other, more well-studied groups (e.g., picture-wing flies) (Hardy et al. 2001). It is unclear whether this suite of traits represents a secondary transition relative to the ancestral state, because it is not known whether haleakalae flies are the sister clade to all other Hawaiian Drosophila or nested within the radiation. Resolution in the relationships at the base of this lineage will be key in identifying which branches experienced substantial trait diversification, and especially in identifying whether any of these traits demonstrate predictable patterns of coevolution.

Here, we present the first phylogenomic relationships between the major groups of Hawaiian Drosophilidae. We combine 12 new transcriptomes sequenced in this study with recently published genomes for two Hawaiian Drosophila species (Kim et al. 2021), four non-Hawaiian Scaptomyza (Kim et al. 2021), and six outgroup species (Larkin et al. 2021). By increasing the number of genes used to infer relationships, we begin to unpack the evolutionary history in the short internodes at the base of the Hawaiian Drosophila radiation. Following up on the critical study by Baker and DeSalle (1997) 25 years ago, we explore the landscape of treespace and the discordance between species and gene trees using our phylotranscriptomic data set. We then use the results of our analysis as initial constraints on subsequent phylogenetic analyses using a data set of 316 species and 44 genes, compiled using all previous phylogenetic studies of Hawaiian Drosophilidae. Finally, we estimate the age of the radiation, and use this time-calibrated tree to identify branches where shifts in trait evolution likely occurred. Our findings suggest a relationship between major clades that is distinct from both previously hypothesized topologies, and that is well supported by both maximum likelihood and Bayesian analyses. We show that examining a comprehensive landscape of tree and trait space can allow for a more complete understanding of evolutionary dynamics in this remarkable island radiation.

Results

Phylotranscriptomics Suggest a New Phylogeny of Hawaiian Drosophilidae

Using a phylotranscriptomic approach, we recovered a new topology between the major clades of Hawaiian Drosophilidae, distinct from those previously hypothesized (fig. 1 and supplementary fig. S1, Supplementary Material online). This topology was the most likely tree estimated using IQtree (Minh, Schmidt, et al. 2020) and RAxML (Stamatakis 2014), as well as the consensus tree with highest posterior probability estimated using PhyloBayes (Lartillot et al. 2013) (fig. 1 and supplementary figs. S2 and S3, Supplementary Material online). Bootstrap support for all branches was 100 and posterior probability was 1, with the exception of the branch subtending the clade uniting MM+AMC (IQtree ultrafast bootstrap: 66, RAxML bootstrap: 57, PhyloBayes posterior probability: 0.52). We also estimated the phylogeny using a multispecies coalescent model with ASTRAL (Zhang et al. 2018), and recovered the same topology with the exception of the placement of D. primaeva (as the sister taxon to PNA, supplementary fig. S4, Supplementary Material online). Each of these analyses were performed on a supermatrix of 10,949 putatively orthologous genes, aligned and assembled using the agalma pipeline (Dunn et al. 2013) with no filtering based on occupancy (actual gene occupancy was 41.7%). To test the senstivity of our results to missing data, we repeated the IQtree analysis on a data set reduced using an occupancy threshold that ensures representation of 80% of taxa at each gene (1,926 genes), and recovered the same topology as with the full set of genes (supplementary fig. S5, Supplementary Material online).

The most likely tree indicates that the PNA clade, including picture-wing species, is the sister clade to all other Hawaiian Drosophila. Drosophilaprimaeva is found to be the sister taxon to a clade containing non-PNA Hawaiian Drosophila, though this clade received lower support when using the data set reduced by occupancy (supplementary fig. S5, Supplementary Material online, ultrafast bootstrap of 85). A second monotypic lineage, D. adventitia, was not sampled for phylotranscriptomic analyses, but using specific gene markers, we recover this as the sister taxon to a clade including MM+AMC+haleakalae (see section on expanded phylogenetic analysis below). This latter clade was recovered in previous phylogenetic analyses (O’Grady et al. 2011; Magnacca and Price 2015). In contrast to those studies, which suggested a monophyletic clade of AMC+haleakalae, we do not recover sufficient support for any particular arrangement of MM, AMC, and haleakalae (ultrafast bootstrap from both the full and reduced occupancy matrix is <95).

We tested the most likely tree emerging from our analysis (fig. 1) against two previously suggested alternative hypotheses (fig. 1) using the Swofford–Olsen–Waddell–Hillis (SOWH) test (Swofford et al. 1996), a parametric bootstrap approach for comparing phylogenetic hypotheses. In both cases, the difference in likelihood between the most likely tree and these alternatives was larger than we would expect by chance (P value for both <0.01, with a sample size of 100). Between figure 1 the difference in log-likelihood was 1,774.1, and between figure 1 was 6,132.1, whereas the null distribution according to the SOWH test had no differences greater than 15 for either comparison. Taken together, our results suggest a new phylogeny for Hawaiian Drosophilidae relationships wherein MM, AMC, and haleakalae represent a monophyletic group, and the PNA clade, rather than either the haleakalae clade or D. primaeva, is the sister clade to all others (fig. 1).

Identifying Hotspots of Gene and Site Concordance in Treespace

We analyzed the strength of phylogenetic concordance in our phylotranscriptomic data set by estimating the gene and site concordance factors for each branch in our tree. Gene concordance factors (gCF) are calculated as the proportion of informative gene trees that contain a given branch between taxa, and can range from 0 to 100 (Baum 2007; Minh, Hahn, et al. 2020). Site concordance factors (sCF) are calculated as the average proportion of informative sites that support a given branch between taxa. Because one site can only support one of three arrangements for a quartet of taxa, sCF typically ranges from ∼33.3 to 100, with 33.3 representing our null expectation based on chance (Minh, Hahn, et al. 2020). We found that for many branches in our tree both gCF and sCF are high, indicating these relationships are supported by a majority of genes and sites in our data set. For example, the branch uniting Hawaiian Drosophila has a gCF of 91.2, and sCF of 72.1 (fig. 1). However, for the branches subtending most relationships between the major clades of Hawaiian Drosophila, gCF and sCF are low. For example, the branch uniting D. primaeva+MM+AMC+haleakalae to the exclusion of PNA has a bootstrap value of 100, but a gCF of 19.3 and sCF of 32.2.

We also tested the extent to which potential error in multiple sequence alignment affected concordance values by filtering out poorly aligned sequence fragments, and repeating the tree inference and concordance analyses. After filtering poorly aligned sequences, we recovered the same topology as figure 1, with the exception of the arrangement of MM, AMC, and haleakalae, here showing a MM+haleakalae as monophyletic (IQtree ultrafast bootstrap: 97, supplementary fig. S6, Supplementary Material online). Concordance factors between analyses on filtered and nonfiltered data were nearly equivalent (e.g., the branch separating PNA from the other Hawaiian Drosophila received a gCF of 19.0 and sCF of 31.6 after filtering, and a gCF of 19.3 and sCF of 32.2 using all data). Furthermore, using a series of stringency thresholds to filter the data, we observed no pattern of increasing or decreasing concordance factors across branches (supplementary fig. S7, Supplementary Material online). These results suggest that values of discordance in this phylogeny are not artificially inflated due to technical errors from the alignment step.

We interpret the measures of discordance as reflecting real variation in the phylogenetic signal of different genes and sites, which is not unexpected for a radiation such as this with short internodes subtending major clades (Minh, Hahn, et al. 2020). Furthermore, the presence of discordance does not mean that there is little that can be said about the relationships between these groups. In contrast, by unpacking this discordance, we can begin to qualitatively describe the amount and distribution of phylogenetic signal for multiple alternative, plausible bipartitions.

To this end, we first visualized hotspots of concordance across treespace (fig. 2). We created all 105 topological combinations of the possible arrangements between major clades, and then re-estimated gCF and sCF for each. Visualizing the mean values for gCF and sCF plotted in treespace shows that the most likely tree, as estimated with IQtree, is not the tree with the highest mean gCF and sCF, but it is near a hotspot of alternative arrangements for which both of these values are high (fig. 2, treespace, most likely tree indicated by dark red outline). In contrast to the most likely topology, the trees with the top three mean gCF values and two of the three trees with the top mean sCF values unite D. primaeva+PNA to the exclusion of other Hawaiian Drosophila. Variation between these top trees largely depends on the placement of haleakalae relative to other clades (fig. 2, top gCF and sCF trees).

Fig. 2.

The landscape of treespace shows hotpots of concordance among genes and sites. The landscape of treespace for all possible topological combinations of the five clades of Hawaiian Drosophila studied here: PNA, D. primaeva, haleakalae, MM, and AMC. Individual points represent different arrangements of the five clades, labeled randomly with two-letter IDs from aa through ea. The distance between points indicates tree similarity (calculated from Robinson–Foulds distances). The size of points represents the mean gene concordance factor (gCF) across relevant branches, and the color represents the mean site concordance factor (sCF; purple, low; yellow, high). The point outlined in red (tree dm) indicates the best topology found with IQtree, RAxML, and PhyloBayes, which is distinct from the top trees according to mean gCF (av, ah, and ap) or mean sCF (av, az, and cf). Concordance measurements for all topologies are available, see Materials and Methods and Data Availability.

Calculating the mean gCF and sCF across branches may not always provide an informative metric, given that some topologies may contain one highly supported branch and others with very low support. Therefore, we also analyzed concordance for all the unique bipartitions across the set of possible topologies (supplementary figs. S8 and S9, Supplementary Material online; see supplementary concordance factor analysis, Supplementary Material online). We found that for gCF, there is clear signal supporting bipartitions that unite D. primaeva+PNA, as well as those that unite MM+AMC+haleakalae (supplementary fig. S8, Supplementary Material online). We found that for sCF, concordance values across bipartitions are more variable, but those that unite PNA+haleakalae show less support than we might expect by chance, whereas those that unite D. primaeva+PNA and AMC+MM show more support (supplementary fig. S9, Supplementary Material online). In addition, between gCF and sCF, we found conflicting signals for bipartitions that define one clade as sister to the rest of Hawaiian Drosophila, with gCF indicating support for PNA (consistent with the most likely topology), and sCF indicating support for haleakalae.

To investigate the source of discordance across genes, we performed a likelihood mapping analysis that assesses the phylogenetic information in each gene (Strimmer and Von Haeseler 1997; Minh, Schmidt, et al. 2020). This analysis calculates the likelihood support for the three possible arrangements of each quartet of taxa in an alignment, and then counts the number of informative quartets that strongly support one arrangement over the other two (Strimmer and Von Haeseler 1997). Here, we performed likelihood mapping on each of the 10,949 genes in our full data set, resulting in 2,182 quartets relevant to the position of the PNA group and 3,075 quartets relevant to the position of the haleakalae group. Our results showed that, in both cases, the vast majority of quartets are uninformative, with no strong likelihood support for any one arrangement (supplementary fig. S10, Supplementary Material online). Although support was for the most part evenly divided among possible arrangements, we observed more quartets that unite haleakalae+MM, to the exclusion of AMC, PNA, and D. primaeva, as well more quartets that unite D. primaeva+PNA to the exclusion of other Hawaiian Drosophilidae, than other arrangements (supplementary fig. S10 and C, Supplementary Material online). We also tested whether genes that support one topology over another at these nodes were enriched for either long or short genes, or fast or slow-evolving genes. Our results showed no correlation between gene tree topology and gene length or evolutionary conservation (supplementary fig. S11, Supplementary Material online), with the exception of genes supporting D. primaeva as the sister to all other Hawaiian Drosophila, which was supported by somewhat slower evolving genes than alternative arrangements (supplementary fig. S11, Supplementary Material online).

In summary, across all analyses, we found consistent evidence for a bipartition that separates PNA and a clade that includes MM and AMC. Although the placement of D. primaeva with MM+AMC+haleakalae received strong support in our maximum likelihood and Bayesian analyses, we observed substantial discordance in this arrangement, and detect signal suggesting a significant amount of shared history between D. primaeva and PNA. Similarly, although the clade uniting MM+AMC+haleakalae received strong bootstrap support, we observed substantial discordance in the placement of haleakalae, and suggest that further resolution in its placement will be possible with additional taxon sampling in that clade.

Calibrating an Expanded Phylogeny to Time

Building on the phylotranscriptomic analyses above, we collected all publicly available genomic and transcriptomic data for species from Hawaiian Drosophila and Scaptomyza. These data were accessioned in nine analyses published since 1997, most of which focused on resolving the phylogenetic relationships within a major clade (Baker and DeSalle 1997; O’Grady and Zilversmit 2004; Bonacum et al. 2005; Lapoint et al. 2011, 2013, 2014; Magnacca and Price 2015; Katoh et al. 2017). The data set we compiled contained 44 genes (six mitochondrial and 38 nuclear) from 316 species (including 271 described and 45 undescribed putative species), with an overall occupancy of 17.3% (supplementary fig. S12, Supplementary Material online). We used this data set to infer the phylogeny with IQtree, constraining the relationships between major clades to conform to the topology shown in figure 1.

The resulting topology is to our knowledge the most species rich phylogenetic tree of the Hawaiian Drosophilidae to date (supplementary fig. S13, Supplementary Material online). Several support values are low (ultrafast bootstrap <95), especially for nodes near the base of the radiation. However, this is not unexpected, given that this phylogeny is estimated primarily from the same data previously analyzed, which recovered alternative relationships at those nodes. Of note are the low support values for the relationships within the MM and haleakalae clades (supplementary fig. S13, Supplementary Material online, polytomies), emphasizing the need for further study in these groups.

We used this expanded genetic data set and topology to estimate the age of the Hawaiian Drosophilidae by calibrating this tree to time using the software package BEAST (Bouckaert et al. 2014). Consistent with recent publications (Obbard et al. 2012; Magnacca and Price 2015; Katoh et al. 2017), our results indicate that the age of the split between Hawaiian Drosophila and Scaptomyza occurred between 20 and 25 million years ago (fig. 3, median root age 22.8 mya). The results shown here were calibrated using updated estimates for the ages at which Hawaiian islands became habitable, based on models of island emergence, growth, and decline via erosion and subsidence (Lim and Marshall 2017) (table 1). Similar results were obtained using the calibration scheme from Russo et al. (2013) that includes a single fossil calibration point for the clade, based on the taxon S. dominicana recovered from dominican amber (median root age 22.9 mya). However, with both schemes, uncertainty around the root age remains substantial (95% highest posterior density confidence interval 17.4–29 mya), and small changes in the calibration times used can lead to substantial differences in this estimate. When calibrating the tree using the same island age estimates as in Magnacca and Price (2015), which are marginally younger (table 1), we estimated the age of Hawaiian Drosophilidae to be ∼15 My old (median root age 15.5 mya). Furthermore, we note that calibrating using primarily vicariance based estimates of time is considered to be imprecise and should be avoided (Kodandaramaiah 2011). Taken together, we consider this estimate of the age of Hawaiian Drosophila, as well as those previously published, to be tentative, and suggest that further data (e.g., new fossil evidence) will be necessary to determine the age of diversification relative to island emergence with greater certainty.

Table 1.

Calibration Points for Dating with BEAST.

Group	Species Included	Magnacca and Price (2015)	Magnacca and Price (2015), Updated Ages	Russo et al. (2013)
Scaptomyza genus	All Scaptomyza species			25.7 (3.0)
planitibia group	D. anomalipes, D. quasianomialipes, D. oahuensis, D. obscuripes, D. hemipeza, D. melanocephala, D. planitibia, D. heteroneura, D. silvestris, D. neoperkinsi, D. neopicta, D. ingens, D. differens, D. substenoptera, D. cyrtoloma, D. hanaulae, D. nigribasis	3 (0.5)	4.135 (0.5)
lanaiensis subgroup	D. lanaiensis, D. hexachaetae, D. digressa, D. moli	3 (0.5)	4.135 (0.5)
picticornis subgroup	D. picticornis, D. setosifrons, D. pilipa			4.7 (0.1)
planitibia subgroup	D. differens, D. hemipeza, D. planitibia, D. silvestris, D. heteroneura			2.8 (0.1)
cyrtoloma subgroup	D. neoperkinsi, D. obscuripes, D. melanocephala, D. cyrtoloma, D. ingens, D. hanaulae, D. oahuensis			2.8 (0.1)
sobrina+orthofascia+ciliaticrus	D. sobrina, D. orthofascia, D. ciliaticrus	1.7 (0.3)	2.55 (0.3)
silvestris+heteroneura	D. silvestris, D. heteroneura	0.5 (0.2)	1.2 (0.2)	0.6 (0.04)

Note.—Values are mean (standard deviation) age in million years for normally distributed time priors.

According to this estimate, we find that the division between major Hawaiian Drosophila clades occurred around 10 mya (fig. 3), prior to the estimated time when the Hawaiian island Kaua’i became habitable (between 6.3 and 6.0 mya; Lim and Marshall 2017). Our results show that the diversification of lineages within MM also occurred around that time, whereas the lineages within the AMC, haleakalae, and grimshawi groups (PNA) all arose within the last 5 million years, around the time Oahu became habitable. We note that the MM groups suffer from lower representation across genes used to calibrate the tree to time (supplementary fig. S12, Supplementary Material online), and suggest that more data may help shed light on differences in the age of this clade relative to others.

Ancestral State Reconstruction of Oviposition and Larval Feeding Ecology

With this time-calibrated tree for 316 species, we have an opportunity to investigate the evolutionary dynamics of trait diversification. By modeling the evolution of the diverse suite of ecological and morphological features across the phylogeny, we can identify which lineages have experienced major shifts in trait evolution. Predicting the number and phylogenetic position of these shifts will in turn be critical for informing future studies on development, life-history, and evolution of these flies. In the following analyses of trait evolution, we used the maximum clade credibility tree from the constrained BEAST analysis described above. Using this tree allows us to maximize the number of taxa for which genetic data are available, painting the most complete picture of ecological and morphological evolution in these flies up to this date. However, due to the fraction of genetic data missing across taxa, it also includes nodes with low bootstrap support (supplementary fig. S13, Supplementary Material online, polytomies). Therefore, for internal lineages for which evolutionary relationships remain uncertain, the position of these evolutionary shifts are subject to change as more genetic data become available and further phylogenetic resolution is achieved.

The Hawaiian Drosophilidae use a wide variety of plant, animal, and fungal species for egg laying and larval feeding (fig. 4) (Heed 1968; Montgomery 1975; Magnacca et al. 2008). The majority of species breed in rotting substrates, with variation in the part of the plant or fungus in question, including rotting bark, leaves, flowers, and fruit. A few species breed on live tissue, and one notable Scaptomyza subgenus, Titanochaeta, have been reared exclusively from spider egg masses (Knab 1914). In 2008, Magnacca et al. (2008) reviewed host plant and substrate records and found that, whereas many species can be considered specialists to species or substrate, host shifting was common and many species occasionally use nonpreferred substrates. The type of oviposition substrate has been suggested as a driver for diversification of the reproductive traits ovariole number and egg size (Kambysellis and Heed 1971; Kambysellis et al. 1995; Sarikaya et al. 2019). However, the previous reconstruction of oviposition substrate by Kambysellis et al. (1995) was performed with a phylogeny that included only three non-PNA species, and was therefore unable to resolve the ancestral oviposition substrate for Hawaiian Drosophila or to identify when evolutionary shifts in substrate outside of PNA were likely to have occurred.

Fig. 4.

Ancestral state reconstruction of oviposition substrate indicates dozens of evolutionary transitions. (A) We used stochastic character mapping to reconstruct the ancestral substrate used for oviposition and larval feeding, and identified dozens of likely transitions in substrate (gray circles). Branch color indicates the ancestral substrate type with highest probability, and tip box indicates extant oviposition substrate. (B) Oviposition substrate category was defined based on rearing records, using the data summarized in Magnacca et al. (2008). Generalist species are defined as those with any two substrates that each comprise of rearing records, or any species with no substrate that comprises of rearing records.

We combined the phylogenetic results presented here with the data summarized in Magnacca et al. (2008), to reconstruct the ancestral oviposition substrate for the Hawaiian Drosophilidae (fig. 5 and supplementary fig. S14, Supplementary Material online). Using stochastic character mapping (Huelsenbeck et al. 2003), we recover the most probable ancestral oviposition substrate for the Hawaiian Drosophila as bark breeding (defined as including rearing records from bark, stems, branches, roots, and fern rachises, see supplementary table S5, Supplementary Material online). We recover a transition from bark to leaf breeding at the base of the AMC clade that has generally persisted throughout the diversification of that group. As previously reported (Magnacca et al. 2008), we find several groups that demonstrate no reported variation in substrate type (e.g., fungus breeding haleakalae, fig. 5).

Fig. 5.

Multiple shifts in evolutionary regimes help explain the diversity of wing length. (A) Using the R package bayou (Uyeda and Harmon 2014), we modeled the evolution of wing length (mm) on the phylogeny and detected several probable shifts in evolutionary regimes (gray circles, larger indicates greater posterior probability that a shift occurred on that branch). Locations of probable shifts include at the base of the MM+AMC+haleakalae clade, subtending the antopocerus group (AMC), and subtending the nudidrosophila (PNA), among others. (B) The distribution of wing lengths across the phylogeny of Hawaiian Drosophila and Scaptomyza.

Over 1,000 stochastic character maps, we recovered an average of 44 transitions in oviposition substrate over the evolutionary history of Hawaiian Drosophilidae. The majority of these changes occurred along branches leading to extant tips, with few transitions at internal nodes (on the summary tree, 8 out of 36 total changes). On average, 70% of transitions were between using a single substrate type as a primary host (“specialist” species) and using multiple types (“generalist” species, defined as using any two substrates that each comprise of all rearing records, or with no substrate that comprises of rearing records; Magnacca et al. 2008). Other transitions were primarily between using rotting bark, leaves, or sap. Pinpointing branches of likely transitions shows that some groups have experienced many more transitions than others, especially MM and non-flower/non-spider egg breeding Scaptomyza. Most generalist species fall in one of these two clades, which also include specialist bark and leaf breeders, among other substrates.

Evolution of Wing, Thorax, and Body Length

Alongside ecological diversification, the Hawaiian Drosophilidae show substantial diversity in adult body size. We used the time-calibrated phylogeny to model the number and timing of major changes in the evolutionary dynamics of size across the phylogeny. First, we digitized 795 records from 26 publications (Grimshaw and Speiser 1901; Grimshaw 1902; Frederick 1914; Bryan 1934, 1938; Wirth 1952; Hackman 1959; Hardy 1965, 1966, 1969, 1977; Hardy and Kaneshiro 1968, 1969, 1971, 1972, 1975, 1979; Kambysellis and Heed 1971; Hardy et al. 2001; O’Grady et al. 2001, 2003; Starmer et al. 2003; Magnacca and O’Grady 2008, 2009; Craddock et al. 2016; Sarikaya et al. 2019), including descriptions of body, wing, and thorax length across 552 species. Then we mapped these traits onto our phylogenetic results, and used the R package bayou (Uyeda and Harmon 2014) to identify branches that represent probable shifts in trait diversification. This package uses Ornstein–Uhlenbeck (OU) models to describe shifts in evolutionary regimes, defined as lineages that share an OU optimum trait value.

In the case of wing length, we find evidence for several highly supported regime shifts in the evolutionary history of Hawaiian Drosophilidae (fig. 5). Some of these are independent shifts on branches subtending groups with larger wings than their nearby relatives, including flies in the antopocerus group (AMC) and in the Engiscaptomyza+Grimshawomyia subgenera (Scaptomyza). Others are independent shifts on branches subtending lineages with smaller wings than nearby relatives such as the nanella+ischnotrix (MM) and the nudidrosophila subgroups (PNA). This suggests that the evolutionary history of Hawaiian Drosophila has included multiple convergent transitions to both larger and smaller wings. In the case of nudidrosophila (PNA), we note that the topology recovered in the summary tree dividing this subgroup into two lineages has very low bootstrap support (supplementary fig. S13, Supplementary Material online, polytomies), and we suggest that the two shifts to smaller wings recovered within PNA may represent a single shift if this group is indeed monophyletic.

We found similar results when considering thorax length (supplementary fig. S15, Supplementary Material online) and body length (supplementary fig. S16, Supplementary Material online). In the case of the former, we find shifts at the base of antopocerus, and nudidrosophila, consistent with the shifts recovered for wing size. In the case of body size, the most probable shifts are located at the base of the Hawaiian Drosophila and the Engiscaptomyza+Grimshawomyia subgenera. However, for body length, no regime shifts received substantially more support than others, despite running bayou for an extra million generations and achieving a final effective size for log-likelihood of 401.9.

Evidence for Convergent Evolution of Ovariole Number and Egg Size

We also performed these analyses on reproductive traits, including egg size, egg shape (aspect ratio, calculated as egg length/width), and the number of egg-producing compartments in the ovary (ovarioles). These traits have been the subject of several life-history studies regarding the hypothesized trade-offs between offspring size and number, and its relationship to ecology (Kambysellis and Heed 1971; Starmer et al. 2003; Sarikaya et al. 2019). Considering egg shape, we find evidence for a shift at the base of the PNA clade, which have proportionally longer eggs than their relatives (supplementary fig. S17, Supplementary Material online). In the case of egg volume, we find evidence for independent shifts on branches subtending flies with large eggs (antopocerus [AMC] and the Engiscaptomyza+Grimshawomyia, supplementary fig. S17, Supplementary Material online). In the case of ovariole number, we find shifts at the base of the MM+AMC+haleakalae clades, which have on an average fewer ovarioles than the other Hawaiian Drosophila (D. primaeva and PNA, supplementary fig. S18, Supplementary Material online).

Work by Kambysellis and Heed (1971) suggested that Hawaiian Drosophilidae species can be grouped into four reproductive categories based on suites of ovarian and egg traits. Subsequent publications (Kambysellis et al. 1995), including work by ourselves (Sarikaya et al. 2019), showed that these categories largely map to differences in oviposition substrate. Given the evidence that ovary and egg traits may be evolving together, we analyzed them with the R package SURFACE (Ingram and Mahler 2013), which uses OU models to analyze regime shifts in multiple traits at once, and allows for distant taxa to share regimes via convergent evolution. The best fitting model indicates four regimes (supplementary fig. S19, Supplementary Material online), two of which correspond to categories defined by Kambysellis and Heed (1971): 1) very large eggs and low ovariole number in S. undulata (Grimshawomyia) and S. nasalis (Engiscaptomyza, group I in their publication); 2) large eggs with moderate to large bodies and moderate ovariole number in antopocerus (AMC) and also in S. crassifemur (Engiscaptomyza) and S. ampliloba (Engiscaptomyza, group II in their publication). The remaining two regimes redistribute species that fall into groups IIIa and IIIb in Kambysellis and Heed (1971) into groups that have 3) small eggs, moderate to large bodies, and high ovariole number in PNA flies, D. primaeva and D. comatifemora (MM); 4) flies with small eggs, small to moderate bodies, and moderate ovariole number, in the remaining MM+AMC flies along with D. preapicula (PNA) and S. devexa (Elmomyza). As predicted by Kambysellis et al. (1995) and ourselves (Sarikaya et al. 2019), these final two regimes are largely divided between bark breeding flies (4) and leaf breeding flies (3).

Discussion

The landscape of treespace, representing support for all the possible topologies given the data, is often hidden from our view (Sanderson et al. 2011; St. John 2017). This is especially true as the size of data sets grows, making it more laborious to traverse treespace landscapes. Approaches such as visualizing the posterior distribution of parameters in a Bayesian analysis, or alternative hypotheses testing (e.g., SOWH test in a maximum likelihood framework), can provide a sense of how support for one result compares with others. But given that a complete exploration of treespace is typically not available, we often do not know whether the support landscape in treespace is generally flat, rugged, or highly structured.

Model clades for phylogenetics such as the Hawaiian Drosophilidae, however, offer an opportunity to explore these methods using real-world data. In the case of the landscape of treespace, especially in the context of discordance of gene trees and species trees, these flies have a long history as one such model clade. Here, we provide a comprehensive snapshot of treespace for this island radiation. We find that, in this case, the landscape of support is largely defined by one hotspot in both gene and site concordance. This hotspot divides the major clades of Hawaiian Drosophila into two main lineages, the picture wing flies and their allies (PNA) on one side, and the modified-mouthparts (MM) and antopocerus, modified-tarsus, ciliated-tarsus (AMC) flies on the other. We consider this division to be strongly indicated given the data, and we note that this is in line with other recent phylogenetic results (supplementary fig. S1, Supplementary Material online) (O’Grady et al. 2011; Magnacca and Price 2015).

Within this hotspot of support, several alternative topologies that differ in the placement of smaller clades (D. primaeva and haleakalae) have an equivocal amount of support across genes and sites. We suggest that much of this discordance represents the results of evolutionary processes that took place on the short internodes at the base of the radiation. Despite this local discordance, the outcome of all phylogenetic software tested here indicates strong support for a single topology (fig. 1). With this information, we consider that tree, with PNA as the sister clade to the rest of Hawaiian Drosophila, and haleakalae as the sister clade to MM+AMC, to be a plausible new hypothesis for the evolution of these flies. We suggest that additional taxonomic sampling in the haleakalae will be valuable in gaining a fine-scale view of the landscape of support within this hotspot.

This new hypothesis for the relationship between major groups has several implications for our understanding of ecological and morphological evolution. Some previous studies have focused on defining one group as “basal” to others (e.g., haleakalae, MM, or D. primaeva) (Kambysellis et al. 1995; O’Grady et al. 2011). However, our results show that an alternative interpretation would serve us better. We find that the PNA clade (including picture-wing flies) is the sister clade to all others, and we note that for at least one trait (bark breeding), most PNA flies appear to display the same state as the last common ancestor of Hawaiian Drosophila. The relationship between this group to haleakalae and others suggests the possibility of a secondary loss of complex courtship behavior in the latter (O’Grady et al. 2011). We note that the overall pattern in the group has been one of many transitions to and from the ancestral state, including in ecology, size, and allometry.

Our results on wing, body, and egg size evolution show that Hawaiian Drosophilidae have experienced multiple, independent shifts to both larger and smaller sizes. These repeated changes present an opportunity to test the predictability of evolution by analyzing whether repeated changes in size are coincident with changes in other features, including ecology, development, and whether these repeated trait changes share the same genetic regulatory basis. The findings of this study on ovariole number and egg size evolution are consistent with what has previously been shown (Kambysellis et al. 1995; Sarikaya et al. 2019), indicating that evolutionary changes in the larval ecology correspond to changes in reproductive trait evolution. However, our findings here show that larval feeding substrate does not explain all the dynamics of trait diversification in Hawaiian Drosophila. For example, the antopocerus group (AMC) shares the same oviposition substrate as most other AMC flies, yet we find evidence that several important shifts in thorax, wing, and egg size evolution all occurred on the branch subtending its diversification.

Previous authors have commented on the potential of the Hawaiian Drosophila as a model clade for the study of the evolution of development (Edwards et al. 2007; O’Grady and DeSalle 2018), given its close relationship to genetic model species like D. melanogaster. Progress in this effort has not always been straightforward, however, given their longer generation times and specific host plant requirements to induce oviposition in the lab (O’Grady and DeSalle 2018). We propose that advances in evo-devo study of the Hawaiian Drosophilidae will be added by leveraging evolutionary methods to formulate and test developmental hypotheses. For example, we can use phylogenetic comparative methods to statistically detect signatures of convergent evolution and to identify changes in patterns of allometric growth (Stevenson et al. 1995; Sarikaya et al. 2019). Going forward, such methods will be essential in providing testable hypotheses regarding the relationship of developmental data to ecological and morphological parameters. The results of these analyses will provide valuable complementary studies to the developmental literature generated using laboratory-amenable model drosophilids, and shed light on the genetic basis of this remarkable island radiation.

Materials and Methods

Field Collection and RNA Extraction

Field Collection

Specimens used for transcriptome sampling were caught on the Hawaiian islands between May of 2016 and May of 2017. Specimens were caught using a combination of net sweeping and fermented banana–mushroom baits in various field sites on the Hawaiian islands of Kaua’i and Hawai’i (see supplementary table S1, Supplementary Material online, for locality data). Field collections were performed under permits issued by the following: Hawai’i Department of Land and Natural Resources, Hawai’i Island Forest Reserves, Kaua’i Island Forest Reserves, Koke’e State Park, and Hawai’i Volcanoes National Park. Adult flies were maintained in the field on vials with a sugar-based media and kept at cool temperatures. They were transported alive back to Cambridge, MA where they were maintained on standard Drosophila media at 18 °C. Samples were processed for RNA extraction between 5 and 31 days after collecting them live in the field (average 9.8 days, see supplementary table S1, Supplementary Material online). One species, Scaptomyza varia, was caught in the field before the adult stage by sampling rotting Clermontia sp. flowers (the oviposition substrate). For this species, male and female adult flies emerged in the lab, and were kept together until sampled for RNA extraction.

Species Identification

Species were identified using dichotomous keys (Hackman 1959; Hardy 1965, 1977; O’Grady et al. 2003; Magnacca and Price 2012) when possible. Many keys for Hawaiian Drosophilidae are written focusing on adult male-specific characters (e.g., sexually dimorphic features or male genitalia) (Hackman 1959). Therefore, for species where females could not be unambiguously identified, we verified their identity using DNA barcoding. When males were caught from the same location, we identified males to species using dichotomous keys and matched their barcode sequences to females included in our study. When males were not available, we matched barcodes from collected females to sequences previously uploaded to NCBI (O’Grady et al. 2011; Lapoint et al. 2014; Katoh et al. 2017).

The following dichotomous keys were used to identify species: for picture-wing males and females, Magnacca and Price (2012); for antopocerus males, Hardy (1977); for Scaptomyza, Hackman (1959); for species in the mimica subgroup of MM, O’Grady et al. (2003); for other miscellaneous species, Hardy (1965).

For DNA barcoding, DNA was extracted from one or two legs from male specimens using the Qiagen DNeasy blood and tissue extraction kit, or from the DNA of females isolated during RNA extraction (see below). We amplified and sequenced the cytochrome oxidase I (COI), II (COII), and 16S rRNA genes using the primers and protocols described in Sarikaya et al. (2019).

For barcode matching, we aligned sequences using MAFFT, version v7.475 (Katoh and Standley 2013), and assembled gene trees using RAxML, version 8.2.9 (Stamatakis 2014). Definitive matches were considered when sequences for females formed a monophyletic clade with reference males or reference sequences from NCBI (supplementary table S2, Supplementary Material online). Sequence files and gene trees are available at the GitHub repository http://github.com/shchurch/hawaiian_drosophilidae_phylogeny_2021, under analysis/data/DNA_barcoding.

Female D. primaeva, D. macrothrix, D. sproati, and D. picticornis could be identified unambiguously using dichotomous keys. Female D. atroscutellata, D. nanella, D. mimica, D. tanythrix, S. cyrtandrae, S. varipicta, and S. varia were identified by matching barcodes to reference sequences from NCBI, reference males, or both. For the female haleakalae fly used in this study, no male flies were caught in the same location as these individuals, and no other sequences for haleakalae males on NCBI were an exact match with this species. Given its similar appearance to Drosophila dives, we are referring to it here as Drosophila cf dives, and we await further molecular and taxonomic studies of this group that will resolve its identity. Photos of individual flies used for transcriptome sequencing are shown in supplementary figure S20, Supplementary Material online.

RNA Extraction

RNA was extracted from frozen samples using the standard TRIzol protocol (http://tools.thermofisher.com/content/sfs/manuals/trizol_reagent.pdf, last accessed February 2022). About 1 ml of TRIzol was added to each frozen sample, which was then homogenized using a sterile motorized mortar. The recommended protocol was followed without modifications, using 10 µg of glycogen, and resuspending in 20 µl RNAse-free water–EDTA–SDS solution. DNA for subsequent barcoding was also extracted using the phenol–chloroform phase saved from the RNA extraction.

RNA concentration was checked using a Qubit fluorometer, and integrity was assessed with an Agilent TapeStation 4200. RNA libraries were prepared following the PrepX polyA mRNA Isolation kit and the PrepX RNA-Seq for Illumina Library kit, using the 48 sample protocol on an Apollo 324 liquid handling robot in the Harvard University Bauer Core Facilities. Final library concentration and integrity were again assessed using the Qubit and TapeStation protocols.

The field collecting for this study was accomplished with a target number of individuals per species in mind, based on future sampling objectives for RNA sequencing studies that, as of the time of writing, have not been published. These objectives were to have four wild-caught, mature, apparently healthy females, three of which were to be dissected for tissue-specific RNA sequencing, and one intended as a whole-body reference library. When four individuals were not available, the reference library was assembled by combining the tissue-specific libraries from one of the other individuals. This was the case for the following species: D. sproati, which was dissected and had RNA extracted separately from the head, ovaries, and carcass, with RNA combined prior to library preparation; and S. varia, S. cyrtandrae, and D. cf dives, for which RNA was extracted and libraries prepared for separate tissues, and raw reads were combined after sequencing.

For the other eight species, sufficient individual females were available for the sampling objectives. In these cases, one entire female fly was dissected and photographed to assess whether vitellogenic eggs were present in the ovary, and all tissues were combined in the same tube and used for RNA extraction.

Libraries for transcriptome assembly were sequenced on an Illumina HiSeq 2500, using the standard version 4 protocol, at 125 base pairs of paired-end reads. A table of total read counts for each library can be found in supplementary table S3, Supplementary Material online.

Transcriptome Assembly

Transcriptome assembly was performed using the agalma pipeline, version 2.0.0 (Dunn et al. 2013). For the 12 new transcriptomes presented in this study, reads from separate rounds of sequencing were concatenated and inserted into the agalma catalog. These were combined with seven publicly available outgroup genomes (D. virilis, D. mojavensis, D. pseudoobscura, D. ananassae, D. willistoni, and D. melanogaster; Larkin et al. 2021), two Hawaiian Drosophila genomes (D. grimshawi, Larkin et al. 2021 and D. murphyi, Kim et al. 2021), and four Scaptomyza genomes (S. graminum, S. montana, S. hsui, S. pallida; Kim et al. 2021). For the non-Hawaiian Drosophila and D. grimshawi genomes, the longest isoform per gene was selected using the gene header. For the four Scaptomyza genomes and D. murphyi genomes, single-copy orthologs were filtered using BUSCO version 4.1.4 (Seppey et al. 2019) against the Diptera obd10 gene set (over 98% of genes were retained as single-copy orthologs). See supplementary table S4, Supplementary Material online, for genome information.

Using the agalma pipeline, the quality score of each library was assessed, and transcriptomes were assembled using the standard parameters. The publicly available genomes were translated and annotated, and the homology of assembled products was inferred using the all-by-all blast component of the homologize step in the agalma pipeline, using nucleotide data and a general time reversible (GTR)+Gamma model to infer gene trees. The agalma version 2.0.0 pipeline also performs a step to reduce the effects of transcript misassignment using a phylogenetically informed approach (treeinform) (Guang et al. 2021). Gene orthology was inferred according to the topology of gene trees estimated with RAxML, orthologs were aligned and trimmed using MAFFT (Katoh and Standley 2013) and Gblocks (Castresana 2000) respectively, and a supermatrix of aligned orthologous sequences was exported.

The final supermatrix output from agalma consisted of 10,949 putatively orthologous genes and 12,758,237 sites. For the primary analyses performed in this manuscript, this supermatrix was not filtered by occupancy, and the actual gene occupancy was 41.9% across the 24 species present in this study. We also created a supermatrix filtered to a target occupancy of 80%, which consisted of 1,926 genes and 1,943,000 sites, which we used to re-estimate the maximum likelihood phylogeny.

All commands used to run the agalma pipeline, and all output report files, are available at the GitHub repository http://github.com/shchurch/hawaiian_drosophilidae_phylogeny_2021, under analysis/phylotranscriptomics.

Phylotranscriptomics and Concordance Factors

We estimated the maximum likelihood phylogeny using IQtree, version 2.1.1 (Minh, Schmidt, et al. 2020). We ran IQtree on a data set partitioned by transcripts, and using the default Model Finder (Kalyaanamoorthy et al. 2017) per partition (Chernomor et al. 2016). For this analysis, partitions containing no informative sites were excluded. We estimated 1,000 ultrafast bootstraps (Hoang et al. 2018) for each node. We also used IQtree to estimate the gene and site concordance factors, as described in Minh, Hahn, et al. (2020). We ran this analysis first on a concatenated data set output by the agalma pipeline command supermatrix with no gene occupancy threshold (returning all aligned transcripts), and then repeated it on a matrix with an 80% occupancy threshold. All subsequent phylogenetic analyses were performed on the full data set.

We also estimated the maximum likelihood phylogeny using the species tree step of the agalma pipeline, which itself runs RAxML, version 8.2.9 (Stamatakis 2014). We used the default parameters for RAxML as called within the agalma phylotranscriptomic pipeline (model GTR+Gamma, 1,000 bootstraps).

We compared the most likely tree against two alternative hypotheses (fig. 1) using the Swofford–Olsen–Waddell–Hillis (SOWH) test (Swofford et al. 1996), as implemented in sowhat, version 1.0 (Church et al. 2015). We ran both comparisons using a GTR+Gamma model, unpartitioned data file, and 100 simulated data sets.

We estimated the phylogeny using the Bayesian software PhyloBayes, mpi version 1.7a (Lartillot et al. 2013). We ran PhyloBayes using a CAT-GTR model for nucleotide data, without partitions, on the full set of transcripts exported from agalma. Phylobayes was run for 1,400–1,900 generations, and convergence was assessed as the maximum difference between two chains. The initial two chains did not show signs of convergence after 1,000 generations (maximum difference was 1), so two additional chains were initialized. These reached convergence with one of the initial chains at 450 generations (maximum difference was 0). The divergence between these three chains and the fourth resulted from differences in the relationships between the MM, AMC, and haleakalae clades. The consensus tree was estimated using all four chains and burn-in of 100 generations, taking every tree (maximum difference was 1).

We estimated the phylogeny using the multispecies coalescent with the software ASTRAL, version 5.7.7 (Zhang et al. 2018). For this analysis, we input the gene trees as inferred by IQtree, using the methods described above.

To further explore the concordance of data across possible topologies in treespace, we wrote a custom python script to create all 105 combinations of possible topologies for the five clades in question, with the root between these clades set at the split between Hawaiian Drosophila and Scaptomyza. We used each of these trees as the constraint for a concordance factor analysis, using the same approach as described above for the most likely tree. We visualized treespace by plotting each tree according to Robinsoun–Foulds distance using the R package treespace, version 1.1.4 (Jombart et al. 2017). We then mapped on this space the mean concordance factors for each topology (calculated as the mean site and gene concordance on branches, excluding those shared between all topologies).

We performed likelihood mapping analyses on each partition from the full supermatrix, using the lmap command in IQtree2 and using all possible quartets. For each partition, we performed two analyses, first analyzing the arrangement of the PNA clade relative to Scaptomyza, D. primaeva, and a clade uniting MM+AMC+haleakalea, and second, analyzing the arrangement of the haleakalae clade relative to MM, AMC, and a clade uniting PNA and D. primaeva. We summarized all likelihood mapping results by counting the total number of quartets across all partitions that fell into each area of a three- and seven-way partitioned likelihood map (Strimmer and Von Haeseler 1997).

We calculated a conservation score for each gene in our supermatrix, using the software trimAL, version v1.4.rev22, as the accumulated similarity score divided by the number of residues in each alignment, following the protocol of Cunha et al. (2021).

We tested the robustness of results to potential errors in multiple sequence alignment by filtering out poorly aligned sequence fragments and repeating the inference and concordance analyses. Poorly aligned fragments were detected and removed using the software spruceup, version 2020.2.19 (Borowiec 2019), with default parameters and cutoffs at 0.95, 0.97, and 0.99 with which outlier fragments were removed.

All commands used to execute the concordance factor analysis are included in the GitHub repository http://github.com/shchurch/hawaiian_drosophilidae_phylogeny_2021 in the rmarkdown file for the supplement of this manuscript, as well as the folder analysis/phylotranscriptomics/concordance-factor.

Estimating an Expanded Phylogeny

We used the phylotranscriptomic results above, combined with previously published genetic data for Hawaiian Drosophilidae, to estimate an expanded phylogeny. First we gathered the accession numbers from all previously published studies of Hawaiian Drosophila and Scaptomyza genetics (Baker and DeSalle 1997; O’Grady and Zilversmit 2004; Bonacum et al. 2005; Lapoint et al. 2011, 2013, 2014; Magnacca and Price 2015; Katoh et al. 2017). Nucleotide data for each accession number were downloaded from NCBI in March of 2019. We made no manual alterations to these sequences, with the following exceptions: We replaced all nonnucleotide sites (e.g., “N,” “R”) with missing data (?); we removed two sequences (U94256.1—D. disjuncta and U94262.1—S. albovittata) from the 16S data set that did not align to other sequences; we manually removed a portion of the COI data set that did not align; we corrected spelling for S. albovittata; and we updated the taxonomic name of D. crassifemur to S. crassifemur. Original and modified sequences are provided at the GitHub repository http://github.com/shchurch/hawaiian_drosophilidae_phylogeny_2021 under analysis/time-calibrated_phylogenetics/downloaded_sequences.

We then aligned these sequences using MAFFT, version 7.457 (Katoh and Standley 2013), –auto option. We visualized alignments, and for gene 16S, we repeated the alignment using the –adjustdirectionaccurately option. We removed all information from the headers with the exception of the species name, and then selected the sequence per species with the fewest gaps. We concatenated sequences using phyutility version 2.2.6 (Smith and Dunn 2008).

Using these concatenated sequences, we estimated a phylogeny for 316 species, including 271 described species and 45 that are undescribed but for which genetic vouchers were available. This tree was estimated using IQtree (Minh, Schmidt, et al. 2020) with the topology constrained using the most likely phylotranscriptomic tree. This constraint tree included only taxa overlapping between the phylotranscriptomic and concatenated data sets, with one exception: D. iki was substituted for D. cf dives, given that this unidentified species was the only representative from the haleakalae clade present in the phylotranscriptomic analysis. No partition model was used for this analysis. We ran IQtree with default parameters, and we estimated 1,000 ultrafast bootstrap support values as well as 1,000 SH-like likelihood ratio tests.

Visualizing the results showed that all major clades (AMC, PNA, MM, haleakalae, and Scaptomyza) were recovered as monophyletic, with the exception of the placement of D. konaensis, a member of the hirtitibia subgroup that was recovered as the sister taxon to the AMC clade. We investigated the source of this discrepancy by analyzing the individual gene trees that had representation for this species (COI, COII, and 16S, tree estimated using IQtree as described above, tree files available at analysis/time-calibrated_phylogenetics/iqtree/iqtree_investigations). These gene trees showed that D. konaensis sequences had variable affinity to unlikely relatives, including Scaptomyza and modified-mouthpart. We considered this to be an artifact of a possible error in accession sequence, and so we removed D. konaensis from downloaded sequences and repeated the alignment and tree estimation steps.

All commands used to download and align sequences as well as estimate the phylogeny, along with all input and output files, are available at the GitHub repository http://github.com/shchurch/hawaiian_drosophilidae_phylogeny_2021, under analysis/time-calibrated_phylogenetics/.

Calibrating the Phylogeny to Time

This expanded phylogeny was calibrated to time using BEAST, version 2.6.3 (Bouckaert et al. 2014). This tree search was performed using the following parameters and priors, set using BEAUTi2 (Bilderbeek and Etienne 2018): a relaxed log-normal clock model, GTR site model with four gamma categories, and a Yule process branching model. We tested three calibration schemes based on previously published analyses: 1) Four normally distributed node priors, used by Magnacca and Price (2015), based on the apparent progression rule seen in the island distribution of these species; 2) the same four node priors, but with island ages adjusted to correspond to recently updated estimates for the age at which islands became habitable (Lim and Marshall 2017), which are based on models that describe the volcanic growth and decay of each Hawaiian island as it has passed over the tectonic hotspot; and 3) four node priors based on progression rule island ages and one node prior based on a single fossil specimen in dominican amber (S. dominicana), used by Russo et al. (2013). Calibration times are shown in table 1.

For all BEAST analysis, the most likely topology from the expanded IQtree search was used to create a starting tree, rooted at the split between Scaptomyza and Drosophila and with branch lengths removed. This topology was fixed throughout the analysis by setting tree topology operator weights to zero.

The BEAST analyses using only island-age based calibrations were run between 20 and 25 million generations. The maximum clade credibility tree was determined using TreeAnnotator (Helfrich et al. 2018), with a burn-in of 10%, selected by visualizing in Tracer, version 1.7.1 (Rambaut et al. 2018). The analysis run using both island-age and the single fossil node calibration was run 50 million generations, and summarized with a burn-in of 5%. For all analyses, the effective size for the posterior was >100 (Magnacca calibration = 581.3, Magnacca updated island ages = 453.7, Russo calibration island ages = 168.7), though for tree height the effective size for the older island ages and fossil-calibrated analyses were both below 100 (Magnacca calibration = 137.4, Magnacca calibration, updated island ages = 92.1, Russo calibration = 87.7).

Estimating Ecological and Morphological Evolutionary Transitions

For ecological data on oviposition and larval feeding substrate, we used the rearing records summarized in Magnacca et al. (2008), Appendix I. Following the method of Magnacca et al., we grouped oviposition substrates into several general categories, listed in supplementary table S5, Supplementary Material online. We also followed the definition from Magnacca and colleagues of nonmonophagous (here referred to as generalist) as any species for which no single substrate type comprised more than of rearing records, or for which any two substrates each comprised more than . We note that D. comatifemora was listed as a bark breeder in Sarikaya et al. (2019), but no rearing records for this species are present in Magnacca et al. (2008, 2009); we therefore listed it here as “breeding habits unknown.”

We reconstructed the ancestral state for general oviposition substrate type using the R package phytools, version 0.7-70 (Revell 2012) on the maximum clade credibility tree from the constrained BEAST analyses. We performed 1,000 simulations of stochastic character mapping using the make.simmap function (with a maximum likelihood method for estimating the transition model), and then summarized the ancestral state at each node as the oviposition substrate with the highest posterior probability. We used this summary tree to identify branches with likely transitions between oviposition substrates.

For morphological data on wing, body, and thorax length, we digitized data from 26 publications (Grimshaw and Speiser 1901; Grimshaw 1902; Frederick 1914; Bryan 1934, 1938; Wirth 1952; Hackman 1959; Hardy 1965, 1966, 1969, 1977; Hardy and Kaneshiro 1968, 1969, 1971, 1972, 1975, 1979; Kambysellis and Heed 1971; Hardy et al. 2001; O’Grady et al. 2001, 2003; Starmer et al. 2003; Magnacca and O’Grady 2008, 2009; Craddock et al. 2016; Sarikaya et al. 2019). For data on ovariole number, egg width, and egg length, we digitized data from three publications (Kambysellis and Heed 1971; Starmer et al. 2003; Sarikaya et al. 2019). We made the following modifications to morphological data: In the data from the GitHub repository associated with the study by Sarikaya et al. (2019), egg size was measured using the radius rather than the diameter; therefore for consistency across studies, we multiplied the reported egg measurements by two. We also excluded data on the egg size of D. incognita from the same publication (Sarikaya et al. 2019) which had two measurements that showed significantly more variation than other measurements (∼150% discrepancy in egg length). We excluded data on wing and thorax length from the manuscript by O’Grady et al. (2003) revising the mimica subgroup, for which all data were inconsistent with measurements from conspecific species in other studies, possibly due to differences in the measurement method.

We identified shifts in evolutionary regimes using the R packages bayou, version 2.2.0 (Uyeda and Harmon 2014) and SURFACE, version 0.5 (Ingram and Mahler 2013) on the maximum clade credibility tree from the constrained BEAST analyses. For all analyses, trait data were log10 transformed. For species that had multiple records for the same trait across publications, we randomly selected one description (data on intraspecific variation or measurement error were not included in analyses due to inconsistency in the methods used to gather and report these data by the original authors). The bayou analyses were performed using a half-Cauchy distribution for the prior value of alpha and sigma2 (scale set at 0.1), a conditional Poisson distribution for the number of shifts (lambda of 10, max of 50), and a normal distribution for theta values (prior mean and SD set at the mean and SD of the trait data). These analyses were run for 1 million generations, with the exception of body and thorax length, which were run for 2 million generations. A burn-in value was set at 0.3 and convergence was evaluated using effective size of the likelihood and the number of shifts (k), see supplementary table S6, Supplementary Material online. The SURFACE analyses were performed on a combination of egg volume, ovariole number, and body length using default parameters and using a two-step forward–backward process of selecting the number of regimes (Ingram and Mahler 2013).

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.