Phylogeography of the Subgenus Drosophila (Diptera: Drosophilidae): Evolutionary History of Faunal Divergence between the Old and the New Worlds

The current subgenus Drosophila (the traditional immigrans-tripunctata radiation) includes major elements of temperate drosophilid faunas in the northern hemisphere. Despite previous molecular phylogenetic analyses, the phylogeny of the subgenus Drosophila has not fully been resolved: the resulting trees have more or less varied in topology. One possible factor for such ambiguous results is taxon-sampling that has been biased towards New World species in previous studies. In this study, taxon sampling was balanced between Old and New World species, and phylogenetic relationships among 45 ingroup species selected from ten core species groups of the subgenus Drosophila were analyzed using nucleotide sequences of three nuclear and two mitochondrial genes. Based on the resulting phylogenetic tree, ancestral distributions and divergence times were estimated for each clade to test Throckmorton’s hypothesis that there was a primary, early-Oligocene disjunction of tropical faunas and a subsequent mid-Miocene disjunction of temperate faunas between the Old and the New Worlds that occurred in parallel in separate lineages of the Drosophilidae. Our results substantially support Throckmorton’s hypothesis of ancestral migrations via the Bering Land Bridge mainly from the Old to the New World, and subsequent vicariant divergence of descendants between the two Worlds occurred in parallel among different lineages of the subgenus Drosophila. However, our results also indicate that these events took place multiple times over a wider time range than Throckmorton proposed, from the late Oligocene to the Pliocene.


Introduction
The genus Drosophila is the most speciose and intensively studied species assemblage in the family Drosophilidae, yet the genus has long been recognized as paraphyletic since the 1975 report of Throckmorton [1], and this characterization has been confirmed in all subsequent, family-wide phylogenetic studies [2][3][4][5][6]. In the most recent molecular phylogenetic analysis [6], the genus Drosophila was considered paraphyletic with respect to 15 different genera. Together, these lineages constitute the tribe Drosophilini and are divided into two major clades: one comprises the genus Lordiphosa and the subgenus Sophophora of Drosophila, and the other includes 14 genera and four subgenera of Drosophila, including the subgenus Drosophila. The subgenus Drosophila is comprised of two main clades [4][5][6][7][8][9][10][11][12][13][14] which are interspersed with eight genera and two other subgenera of Drosophila [6]. One clade of the subgenus Drosophila had been referred to as the virilis section [15] or the virilis-repleta radiation [1,16], and the other as the quinaria section [15] or the immigrans-tripunctata radiation [1,16]. To resolve this paraphyly, Yassin [6] revised the phylogenetic classification of the subgenus Drosophila: the virilisrepleta clade (the virilis section) was revised as the subgenus Siphlodora, while the immigranstripunctata clade (the quinaria section) was left in the subgenus Drosophila. The revised subgenus Drosophila includes 21 species groups [6], of which 15 have been regarded as members of the immigrans-tripunctata clade [4]. They are widely distributed mainly in the northern hemisphere, with four species groups endemic to the Old World (the ancora, bizonata, histrio and immigrans groups), 14 to the New World (the appendiculata, calloptera, cardini, ecuadoriensis, guarani, guttifera, lutzii, macroptera, pallidipennis, peruensis, rubrifrons, sticta, tripunctata and xanthopallescens groups) and three distributed in both (the funebris, quinaria and testacea groups). Some of them, such as the funebris, guttifera, quinaria and testacea groups, represent major elements of temperate drosophilid faunas in the northern hemisphere.
The revised subgenera Siphlodora and Drosophila show parallel biogeographic patterns, i.e., the disjunctions between the Old and the New Worlds in tropical elements (species groups) and in species of temperate species groups. Throckmorton [1] hypothesized that this overall pattern of faunal disjunction between the Old and the New Worlds had emerged in parallel in five lineages (radiations in his descriptions) of Drosophilidae, i.e., in the Scaptodrosophila radiation, the sophophoran radiation, the virilis-repleta radiation (the revised Siphlodora), the immigrans-tripunctata radiation (the revised Drosophila) and the Hirtodrosophila radiation. However, this hypothesis was not presented in a testable way based on solid evidence of phylogenetic relationships. To remedy this, multiple DNA sequence-based studies have provided estimates of divergence times and biogeographic analyses of Old World and New World groups in some lineages of the Drosophilidae [17][18][19][20].
Phylogenetic relationships within the revised subgenus Drosophila or the immigrans-tripunctata radiation have been analyzed in a number of studies based on molecular data [2,3,5,6,[12][13][14][19][20][21][22][23][24]. However, the resulting trees have varied in topology, especially with respect to the basal branches, and have not supported the monophyly of some species groups. Two recent molecular phylogeographic analyses failed to fully resolve the relationships in the immigrans-tripunctata radiation with high confidence [19,20]. With less resolved phylogenetic hypotheses, estimates of divergence times and ancestral distributions can become equivocal, and inferences for evolutionary history are controversial. Limited taxon sampling can contribute to these ambiguous results, and in many cases such sampling has been biased to New World taxa [18].
We gathered species widely from both the New World and the Old World, selecting 22 species from the Old World, 21 from the New World, one from the Australasian region and one cosmopolitan species in an attempt to reduce the taxon sampling bias. We estimated phylogenetic relationships by analyzing the nucleotide sequences of three nuclear and two mitochondrial genes and estimated ancestral distributions and divergence times in order to test Throckmorton's [1] hypothesis about parallel patterns of Old World-New World divergence in the evolution of Drosophilidae.

Taxon sampling
A total of 45 species were selected as focal ingroup taxa from ten core species groups (bizonata, cardini, funebris, guarani, guttifera, histrio, immigrans, quinaria, testacea and tripunctata) of the revised subgenus Drosophila (the immigrans-tripunctata radiation). As outgroup taxa, D. virilis, D. hydei and Idiomyia grimshawi were selected from the sister clade of the subgenus Drosophila [6,19,20]. DNA was extracted from either live or alcohol-preserved specimens. The specimens of D. innubila, D. recens, D. tenebrosa and D. transversa were kindly provided by Prof. John Jaenike (University of Rochester). The remaining 41 specimens were obtained from the UC San Diego Drosophila Stock Center (DSSC), the Tokyo Metropolitan University (TMU) or field collections (Table 1). We determined the DNA sequences of alcohol dehydrogenase (Adh) for 35 species, glycerol-3-phosphate dehydrogenase (Gpdh) for 36 species, mitochondrial 28S ribosomal RNA (28S) for 45 species, cytochrome c oxidase subunit I (COI) for 45 species, and cytochrome c oxidase subunit II (COII) for 45 species. GenBank accession numbers of these sequences, along with those for known sequences used in the analyses, and taxonomic and biogeographic information are given in Table 1.
PCR thermal cycling conditions were 1 min at 94°C; 35 cycles of 30 s at 94°C, 30 s at annealing temperature and 90 s at 72°C; and 7 min at 72°C. Annealing temperatures were 45°C (COI), 49°C (Adh), 51°C (COII), 55°C (28S), or 60°C (Gpdh). The primers used for the amplifications are listed in Table 2. Almost all of the PCR products were sequenced directly. However, 27 Adh and nine Gpdh sequences were sequenced after cloning into pGEM-T Easy vector (Promega) using E. coli DH5α as the host. Sequences were determined in both directions using a BigDye Terminator Sequencing Kit (Life Technologies) and an ABI 3130 and 3730 Genetic Analyzer, according to the manufacturer's protocols.

Phylogenetic analyses
Nucleotide sequences were aligned with the MUSCLE [28] algorithm implemented in MEGA 5.05 [29] with the default settings: gap opening = -400, gap extend = 0, max iterations = 8, clustering method = UPGMB and minimum diagonal length = 24. The intron sequences of Adh and Gpdh were completely removed. Gaps were produced in the alignment procedure. The aligned sequence of 28S only included 15 gap sites. The resulting alignments were checked by eye. The aligned sequences of the five genes were concatenated by using FASconCAT [30] and the total length of the concatenated alignment was 3,443 bp (S1 Dataset).
We used the maximum-likelihood (ML) and Bayesian inference (BI) methods to construct phylogenetic trees. Partitioning concatenated data into more homogeneous subsets improves parameter estimation, reducing small data subset effects that can reduce signal-to-variance ratios [31]. Therefore, the coding sequences of each gene were partitioned into the 1st+2nd and 3rd codon positions, but the 28S sequences were not partitioned. The optimal substitution models for ML and BI were selected by using jModelTest 2.2 [32,33] on the basis of the corrected Akaike information criterion (AICc) [34]. Partitions and their associated substitution models are shown in Table 3. The ML analysis was performed using RAxML 8.0.26 with RAxMLGUI 1.3.1 [35,36]. To obtain the ML tree, we used the rapid hill-climbing algorithm [37] with the 1,000 slow bootstrap (BP) option. The evolutionary model used for each data partition was GTR+I+G.
MrBayes 3.2.3 [38,39] was used to obtain the BI tree. For BI analyses, a starting tree was randomly selected and four chains were run. In each BI analysis, two independent runs were implemented in parallel, with trees sampled every 1,000 cycles. Runs were stopped after 5,000,000 Markov Chain Monte Carlo (MCMC) cycles. For each run, 5,001 samples were obtained, among which 2,500 early-phase trees were discarded as burn-in. The trace files generated by the MCMC runs were inspected in TRACER 1.6.0 [40] to check whether the number of sampling generations and effective sample sizes (ESS) were large enough for reliable parameter estimates. A consensus of sampled trees was computed using the"sumt" command, and the posterior probability (PP) for each interior branch was obtained to assess the robustness of the inferred relationships.

Estimation of divergence times
Divergence times were estimated using the relaxed molecular clock method implemented in BEAST 2.1.3 [41]. This method permits evolutionary rate variation among lineages. An uncorrelated lognormal rate evolutionary model was assumed for descendant lineages, and MCMC was used to estimate the parameters. The alignment and evolutionary models of nucleotide substitutions were the same as those used in the phylogenetic tree construction. A Yule prior was applied for the speciation process. Four MCMC chains of 50 million generations each were sampled every 1,000th generation. We collected 75% post burn-in samples and used TRACER to check whether the number of generations and ESS were large enough for reliable parameter estimation. No fossil records that can be used as calibration points exist in the subgenus Drosophila, hence we adopted the following estimates as calibration points from Russo et al. [20]: the estimated crown ages for the revised subgenera Drosophila (34.5 ± SD 12.3 Mya) and Siphlodora (27.1± 4.6 Mya) were introduced as priors with normal distributions in the analysis.

Biogeographical analysis
Distribution ranges of extant species were characterized by six unit biogeographic areas (OR: Oriental; PA: Palearctic; NA: Nearctic; NT: Neotropical; AU: Australasian; and AF: Afrotropical) Table 2. List of primers for the target genes.

Locus Length (bp) Primers
Adh 711 Adh-F: 5'-ATGGCAATCGCTAAGAA-3'  and their combinations, referring to DrosWLD [42]. In order to infer possible ranges of ancestral species, the Bayesian Binary MCMC (BBM) method implemented in RASP 3.0 [43] was used. This method is primarily designed for reconstructing the ancestral state of a given node. In the BBM analysis, the probabilities of ancestral ranges at a node were estimated in terms of the MrBayesgenerated probabilities of the unit biogeographic areas or their combinations. For estimating the possible ranges of ancestral species, 50,001 binary trees generated by BEAST 2.1.3 were used as input trees, with 12,501 discarded as burn-in. The maximum number of possible ancestral ranges at each node was set to six and the model of BBM analysis was as follows: the state frequency was fixed (Jukes-Canter) and the among-site rate variation was equal. For inferring possible ranges of ancestral species, 10 chains were run simultaneously for 50,000 generations of MCMC. The state was sampled every 100 trees and 100 early-phase trees were discarded as burn-in. In Figs 1 and 2, the ingroup species of the subgenus Drosophila were divided into two major clades, the immigrans species group and a clade that included all the remaining species groups. Within the immigrans group, the basal topology was not highly resolved except for the strongly supported monophyly of the nasuta species subgroup including D. albomicans, D. nasuta, D. neonasuta and D. pallidifrons. Although the immigrans subgroup (D. immigrans, D. formosana and D. ruberrima) and the hypocausta subgroup (D. hypocausta, D. neohypocausta and D. siamana) also formed monophyletic groups, the support for these clades was not as high, especially in the ML tree (Fig 1).

Phylogeny
The other clade comprised 35 species belonging to nine species groups (bizonata, cardini, funebris, guarani, guttifera, histrio, quinaria, testacea and tripunctata). Among the other species groups that included multiple representatives, only two, the funebris group and the testacea group, were recovered as highly supported clades. The quinaria group was paraphyletic, but formed a highly supported clade along with the monotypic guttifera group. The clade including the quinaria + guttifera groups was split into two subclades with high confidence; one was  (Figs 1 and 2 and S1 Fig). The two species of the histrio group were placed distantly: D. sternopleuralis was placed as the most basal branch to the clade comprising all the   Table 1. The red typeface numbers indicate node numbers corresponding to those in Table 4. doi:10.1371/journal.pone.0160051.g002 Phylogeography of the Subgenus Drosophila (Diptera: Drosophilidae) remaining 34 species, but D. histrio was assigned as the sister branch to the quinaria + guttifera clade. D. bizonata (the bizonata group) formed a highly supported clade with the testacea group. The three species groups (tripunctata, cardini and guarani) endemic to the New World formed another highly supported clade. Within this clade, however, the cardini group was not monophyletic. One of its members, D. acutilabella, formed a clade with D. tripunctata (the tripunctata group). However, the confidence for the basal relationships among four lineages (nodes 9, 12, 17 and 20 in Fig 2) was not very high.  Table 4). From this ancestor, two lineages diverged in different geographic regions: the MRCA of the immigrans group was most likely distributed in OR (76%) at about 20.7 Mya (node 10), while the MRCA of the other lineage including the bizonata, cardini, funebris, guarani, guttifera, histrio, quinaria, testacea and tripunctata groups in PA (78%) was likely distributed at about 31.0 Mya (node 4). The most recent range expansion was seen in D. nasuta (node 43), which very recently invaded the Neotropical (NT) region from its original range of AF and OR [44]. Dispersals from the Old World to the New World were inferred to have occurred five times in the latter lineage. The first migratory event occurred from the PA (45%) or PA+OR (38%) distribution (node 8) to Nearctic (NA)+NT (39%) or NT (38%) (node 12) in the tripunctata, cardini and guarani group clade. The second event was detected in the ancestor of the funebris group, which migrated from PA (67%) or PA+OR (29%) to the NA (71%) during the period of 26.7 (node 5) to 15.5 (node 17) Mya. The third event was a migration or range expansion in the ancestor of the Nearctic quinaria-group clade including D. guttifera, which occurred at 19.5 (node 11) to 9.6 (node 22) Mya. The fourth was a migration in the ancestor of D. falleni and D. innubila from PA (53%) or PA+OR (40%) to NA (93%) 9.5 (node 23) to 4.3 (node 34) Mya. The fifth was a migration of D. putrida from PA (57%) or PA+OR (32%) to NA (current range) after 8.0 Mya (node 25).

Biogeographic history
In addition, dispersals in the reverse direction from the New World to the Old World were inferred in a few individual species and a small clade. From the above-mentioned data, Nearctic MRCA of the funebris group, the ancestor of the Old World species D. funebris and D. multispina, migrated back to PA (84%) or PA+OR (16%) 26.7 (node 17) to 15.5 (node 24) Mya. Within the Nearctic clade of the quinaria group, there are two species that are not endemic to NA: D. nigromaculata distributed in PA and D. transversa distributed in the Holarctic (PA +NA). The PA distribution of D. nigromaculata would have resulted from the backward migration of this species from the NA ancestor or vicariant speciation of the Holarctic ancestor into this PA species and the NA ancestor of the other species. However, this inference was not supported in the ML tree (Fig 1) or BI tree (S1 Fig): D. nigromaculata was placed basal to the clade containing D. guttifera and the other six species of the quinaria group. D. transversa expanded its range to PA after 3.0 Mya (node 40).

Phylogeny
The phylogeny reconstructed by the present study was largely consistent with those in previous studies. The subgenus Drosophila (the traditional immigrans-tripunctata radiation) was  consistently divided into two clades: one was the immigrans species group and the other comprised the remaining species groups [2,5,13,14,20,23,24]. Although the immigrans group has thus far been regarded as monophyletic, this may be due to insufficient taxon-sampling for this species group. Recent studies have revealed that D. quadrilineata, a member of the quadrilineata subgroup of the immigrans group, was placed outside the clade of the immigrans group proper [4,12,19], and that it is closely related to the genus Samoaia [6,45]. Moreover, Katoh et al. [12] suggested that D. annulipes, another species of the quadrilineata subgroup, is related to the genus Idiomyia (or Hawaiian Drosophila), along with D. maculinotata. Thus, the quadrilineata subgroup is still questionable in terms of its monophyly and phylogenetic position [12]. Another questionable species subgroup in the immigrans group is the curviceps subgroup. When Zhang and Toda [46] established this species subgroup, they suggested that it is closer to the quadrilineata subgroup than to the other subgroups of the immigrans group. Yassin [6] did not assign these two species subgroups to any subgenus in his revised classification of Drosophila. This taxonomic situation is inconvenient: some members (the immigrans, nasuta and hypocausta subgroups) of the immigrans group belong to the subgenus Drosophila, but the others (the quadrilineata and curviceps subgroups) are incertae sedis for the subgenus. To partly solve this problem, Pradhan et al. [47] have made the curviceps subgroup independent, as a species group, from the immigrans group, but left its subgeneric assignment undetermined, i.e., incertae sedis. Thus, more taxonomic revision is needed for the traditional immigrans group and its related taxa. The clade including the species groups other than the immigrans group is here termed the quinaria-tripunctata radiation and is characterized by the less resolved, basal branching of major component clades. Within this radiation, two large, well-supported clades were recognized in agreement with previous studies. One is the so-called tripunctata radiation, the monophyly of which has been repeatedly inferred in the previous studies [13,14,[22][23][24]. Currently, this radiation includes 123 species of six (calloptera, cardini, guarani, pallidipennis, sticta and tripunctata) species groups that diversified mostly in the Neotropical region [22,24,42]. Phylogenetic relationships among these species groups and the monophyly for some of them are still questionable [5,13,14,23,24], and the present study cannot address these issues because of poor taxon-sampling (only 11 species from three species groups). The other large, well-supported clade comprises the quinaria and guttifera species groups. It has often been indicated in previous analyses [4,5,13,[19][20][21]48] that D. guttifera is closely related to the quinaria group. Accordingly, Markow and O'Grady [49] treated D. guttifera as a species of the quinaria group. The present study has confirmed the inference by Perlman et al. [21] that this "quinaria" clade is further divided into two subclades. These subclades can be more or less distinguished from each other by their choice of breeding niches: one consists of fungal breeders such as D. phalerata, D. brachynephros and D. falleni, whereas the other includes decaying plant breeders such as D. nigromaculata, D. deflecta, D. quinaria, D. limbata, D. subpalustris and D. palustris [21,50,51].
However, some phylogenetic inferences drawn from the present study are new or different from those of previous studies. The phylogeny of the histrio species group has been unclear for a long time, although either D. histrio or D. sternopleuralis has been included in some previous studies [3][4][5][6]20,21]. The present study included both species and revealed that the histrio group is not monophyletic. With high confidence, D. sternopleuralis represents the most basal lineage in the quinaria-tripunctata radiation, while D. histrio is the sister taxon to the "quinaria" clade. The phylogenetic position of D. sternopleuralis has been problematic in previous studies [3,4,20]. It should be noted, however, that Russo et al. [20] placed D. sternopleuralis along with D. pruinosa of the Afrotropical loiciana species complex [52] in the sister clade to the immigrans group proper, though not highly supported. Furthermore, although D. sternopleuralis is distributed in warm temperate regions of Japan and Korea, a number of its relatives have been recorded in the tropical Australasian to the subtropical Oriental regions, i.e., Papua New Guinea, Borneo, Sri Lanka, India, Nepal, Myanmar, southern China and Taiwan [43,[53][54][55][56]. Although the phylogenetic position of D. histrio has also been controversial [4][5][6]21], our results strongly support that it is a sister relationship to the "quinaria" clade, which is consistent with the topology inferred by van der Linde and Houle [4].
Only two previous studies [5,6] examined the monophyly of the funebris species group, including D. funebris and D. macrospina, in molecular phylogenetic analyses. However, the results were controversial. The present study suggests that three species, D. multispina in addition to the two above, of the funebris group form a highly supported clade. Although its placement within the quinaria-tripunctata radiation is consistent with most previous studies [2,[3][4][5][6]9,14,[19][20][21][22], its phylogenetic position within the radiation has been uncertain in previous studies.
The testacea and bizonata species groups formed a clade in previous, phylogenetic studies [5,6], and our analyses confirm this, but differ from that of van der Linde et al. [5] with respect to the phylogenetic position of D. bizonata. We placed it as the sister to the testacea group, which is monophyletic. However, van der Linde et al. [5] placed D. bizonata within the testacea group, rendering the latter paraphyletic. The discrepancy may have been due to the difference in gene markers employed for the phylogenetic analyses. In addition, both studies included only one species (D. bizonata) of the bizonata group, of which six more species are known. This issue should be addressed using better taxon-sampling.

Evolutionary history
Throckmorton [1] hypothesized that two major biogeographic patterns resulted in the evolution of the Drosophilidae. The first is the primary disjunction of tropical faunas between the Old and the New Worlds in the early Oligocene, and the second is the disjunction of temperate faunas between the two Worlds in the mid-Miocene. A number of recent studies have tested this hypothesis by molecular phylogenetic analyses in specific lineages [17,18,24] or across large taxa such as the genus Drosophila [19], the subfamily Drosophilinae [23] and the entire family Drosophilidae [20]. However, the time estimates for these disjunctions varied among the studies because of differences in the timescale placed on the phylogenies that were used [57]. Dates from two sources have been used as calibration points to estimate divergence times within the Drosophilidae: one from fossils, and the other from biogeography. The former source is rare in Drosophilidae, with only few Baltic and Dominican amber fossils dated geologically [58,59]. The latter source comes from the phylogeography of Hawaiian drosophilids, and has more often been used to establish calibration points in molecular phylogenetic studies on Drosophilidae [60,61]. However, Obbard et al. [57] pointed out some biases of the Hawaiian-calibrated dates. Secondarily, taxon sampling affects the inferences on evolutionary history through estimations of ancestral distributions and divergence times. Gao et al. [18] pointed out the importance of recognizing genuine sister groups in obtaining reliable estimates for these attributes.
Since neither fossil-based nor phylogeography-based calibration points are available for the species employed in the present study, we applied the divergence time estimates at two nodes (node 2 and node 6) from Russo et al. [20] as calibration points. The study by Russo et al. employed the largest number of calibration points (two fossil-based and six Hawaiian-phylogeographic calibrations) for estimating divergence times, and was based on the largest taxon samplings (358 species) from the Drosophilidae, but it was somewhat biased to New World species for the quinaria-tripunctata radiation. The present study was conducted to compensate for this bias by including more Old World species of this radiation.
Russo et al. [20] inferred that the pattern of drosophilid faunal disjunctions between the Old and the New Worlds was not as simple as Throckmorton [1] depicted. Such disjunctions were regarded as having occurred multiple times during the long geological period from the late Eocene to the Holocene. However, the most important aspect of Throckmorton's [1] hypothesis, i.e., the parallelism of disjunctions on separate lineages of Drosophilidae, was recovered, indicating that such disjunctions may have been caused by vicariant events. The present study, focusing on the revised subgenus Drosophila (the immigrans-tripunctata radiation), has corroborated the general pattern of disjunctions inferred by Russo et al. [20], but has shown that the evolutionary history of this lineage should largely be revised.
The first difference was seen in the estimation of ancestral distribution for the quinaria-tripunctata radiation. Russo et al. [20] estimated that the MRCA of this radiation was most likely distributed in the New World, while our results suggested that it was distributed in the Palearctic region. Certainly, this could be due to the difference in taxon sampling between the two studies. Accordingly, Russo et al. [20] inferred that the first Old World-New World disjunction within the subgenus Drosophila occurred in the late Eocene to produce two tropical lineages, the immigrans species group in the Oriental region and the tripunctata radiation in the Neotropical region. In parallel with this disjunction, another vicariant divergence of tropical/subtropical elements was estimated to have concurred in the late Eocene between the Oriental/ Palearctic genus Lordiphosa and the Neotropical saltans and willistoni species groups of the subgenus Sophophora. However, our results suggest that the ancestor of the quinaria-tripunctata radiation had once entered the temperate Palearctic region in the early to mid-Oligocene. From this ancestor, then, two lineages, the funebris species group and the tripunctata radiation, arose by nearly concurrent migrations of their ancestors from the Palearctic region to the New World in the late Oligocene to the early Miocene. Gao et al. [18] estimated that the primary, tropical disjunctions such as that between Lordiphosa and Neotropical Sophophora occurred in the mid-Eocene when the North Atlantic Land Bridge (NALB) broke. Even if the NALB could serve as a migrating route for organisms such as drosophilid flies with larger dispersal abilities until later ages [20,62], our estimate of the ancestral age of the Neotropical tripunctata radiation was much younger than the time of the NALB break. In addition, this lineage was regarded as derived from an ancestor distributed in the Palearctic and/or the Oriental regions in the late Oligocene, when tropical forests continuously graded into temperate forests only in East Asia [62,63]. This suggests that the ancestor of the tripunctata radiation, along with the ancestor of the funebris group, migrated through the Bering Land Bridge into the Nearctic region. Robe et al. [24] inferred, based on molecular phylogenetic and biogeographic analyses, that after the colonization to the Nearctic region this lineage had expanded its range to the Neotropical region and much diversified there.
The second difference was seen in the estimated time and direction of ancestral migration between the Old and the New Worlds in the "quinaria" clade. Russo et al. [20] estimated that the ancestor of this clade had migrated from the New World to the Old World in the late Oligocene. Our divergence time estimates (Fig 2) suggest that it is likely that the ancestor was distributed in the Oriental and/or the Palearctic regions. From this Old World ancestor, one lineage was derived by migration or range expansion to the Nearctic region in the mid-Miocene. Another species concurrently migrated, but in the opposite direction, from the Nearctic to the Palearctic region and became the ancestor of Old World species in the funebris group. The Old World-New World temperate disjunctions dated to the mid-Miocene age were detected in the lineages of the genus Hirtodrosophila and the Drosophila (Siphlodora) melanica species group as well [20].
The third phase of Old World-New World disjunctions was estimated to have occurred on three lineages in the "quinaria" clade and the testacea species group, as well as in the Drosophila (Sophophora) obscura species group [20], in the late Miocene. In the fourth phase, the subarctic species D. transversa achieved its Holarctic distribution by expanding its range from the Nearctic to the Palearctic region in the Pliocene. In this phase, cool-temperate species of the Drosophila (Siphlodora) robusta and virilis species groups became disjunct between the Palearctic and the Nearctic regions [20].
Thus, ancestral migrations and subsequent vicariant divergences of descendants between the Old and the New Worlds occurred multiple times in parallel among different lineages in the quinaria-tripunctata radiation of the subgenus Drosophila. Such faunal disjunctions between the two Worlds concurred in other drosophilid lineages as well, and shaped a part of the overall phylogeographic pattern in the family Drosophilidae [20]. These faunal vicariant events correspond to the secondary disjunctions of temperate elements in Throckmorton's [1] hypothesis, but would have occurred repeatedly in a wider time range, from the late Oligocene to the Pliocene, than Throckmorton [1] proposed. Such multiple migration episodes between the Old and the New World have also been well documented in terrestrial plants and mammals. Tiffney [64] proposed that floristic exchanges between the two Worlds most likely occurred in five major periods: (1) the pre-Tertiary, (2) the early Eocene, (3) the late Eocene-Oligocene, (4) the Miocene, and (5) the late Tertiary-Quaternary periods. The four phases of Old World-New World disjunctions detected for the quinaria-tripunctata radiation in the present study seem to correspond to the last two of the five flora-exchange periods. Zhanxiang [65] recognized three major waves of dispersals from Eurasia to North America in carnivoran mammals: the first in the early Miocene, the second in the late Miocene and the third in the Pliocene: corresponding to the first, third and fourth phases of faunal disjunctions in the quinaria-tripunctata radiation, respectively. These parallel, spatiotemporal patterns observed in terrestrial organisms of the northern hemisphere strongly suggest the contribution of the Bering Land Bridge as a common factor to exchanges and isolation of temperate elements between the Old and the New Worlds [20,62,65]. Paleo-climatic and -vegetational conditions around Beringia gradually changed with global cooling after the Early Eocene Climatic Optimum [66]. The vegetation around the Beringia was dominated by rich deciduous and some coniferous temperate forests in the late Oligocene to the early Miocene, but changed to the cool-temperate forest with decreased diversity of deciduous angiosperms and conifers thereafter from the late Miocene to the Pliocene and to the taiga/tundra in the Holocene [62,63]. The changing climatic and vegetational conditions around the Beringia should have filtered different elements of the temperate biota in the process of faunal/floral exchanges between the Old and the New Worlds after the early or mid-Miocene age: elements adapted to the warmer climate became disjunct earlier in geological time.
Our results clearly show the importance of wider and denser taxon sampling for the reconstruction of evolutionary histories, although there still remain a number of species groups of the subgenus Drosophila to be included in phylogenetic analyses. Further expansion of taxon sampling across the whole family Drosophilidae should provide more evidence for parallel biogeographic disjunctions in separate lineages, and will help to document the vicariant nature of regional faunal development.
Supporting Information S1 Dataset. The nexus format file of sequences used in the analyses.