What does mitogenomics tell us about the evolutionary history of the Drosophila buzzatii cluster (repleta group)?

The Drosophila repleta group is an array of more than 100 species endemic to the "New World", many of which are cactophilic. The ability to utilize decaying cactus tissues as breeding and feeding sites is a key aspect that allowed the successful diversification of the repleta group in American deserts and arid lands. Within this group, the Drosophila buzzatii cluster is a South American clade of seven closely related species in different stages of divergence, making them a valuable model system for evolutionary research. Substantial effort has been devoted to elucidating the phylogenetic relationships among members of the D. buzzatii cluster, including molecular phylogenetic studies that have generated ambiguous results where different tree topologies have resulted dependent on the kinds of molecular marker used. Even though mitochondrial DNA regions have become useful markers in evolutionary biology and population genetics, none of the more than twenty Drosophila mitogenomes assembled so far includes this cluster. Here, we report the assembly of six complete mitogenomes of five species: D. antonietae, D. borborema, D. buzzatii, two strains of D. koepferae and D. seriema, with the aim of revisiting phylogenetic relationships and divergence times by means of mitogenomic analyses. Our recovered topology using complete mitogenomes supports the hypothesis of monophyly of the D. buzzatii cluster and shows two main clades, one including D. buzzatii and D. koepferae (both strains), and the other containing the remaining species. These results are in agreement with previous reports based on a few mitochondrial and/or nuclear genes, but conflict with the results of a recent large-scale nuclear phylogeny, indicating that nuclear and mitochondrial genomes depict different evolutionary histories.

Introduction

Almost all mitochondrial genomes, the “mitogenome”, can be assembled directly from genome or even transcriptome sequencing datasets [1, 2]. The exponential development of next-generation sequencing (NGS) technologies, together with efficient bioinformatic tools for the analysis of genomic information, has allowed efficient assembly of mitochondrial genomes, giving rise to the emergence of the mitogenomics era [3]. Mitogenomics has been very useful in illuminating phylogenetic relationships at various depths of the Tree of Life, e.g., among early branching of metazoan phyla [4], among crocodilians and their survival at the Cretaceous-Tertiary boundary [5], primates [6], the largest clade of freshwater actynopterigian fishes [7] and Anura, the largest living Amphibian group [8]. Also, mitogenomic approaches have been used to investigate evolutionary relationships in groups of closely related species (e.g. [9]). In animals, the mitochondrial genome has been a popular choice in phylogenetic and phylogeographic studies because of its mode of inheritance, rapid evolution and the fact that it does not recombine [10]. Such physical linkage implies that all regions of mitogenomes are expected to produce the same phylogeny. However, the use of different mitogenome regions or even the complete mitogenome may lead to incongruent results [11], suggesting that mitogenomics sometimes may not reflect the true species history but rather the mitochondrial history [12-16]. Moreover, there is evidence suggesting that mtDNA genes are not strictly neutral markers, casting doubts on its use to infer the past history of taxa [17]. Inconsistencies across markers may result from inaccurate reconstructions or from actual differences between genes and species trees. In fact, most methods do not take into consideration that different genomic regions may have different evolutionary histories, mainly due to the occurrence of incomplete lineage sorting and introgressive hybridization [18-20].

Over the last century, the Drosophila genus has been extensively studied because of the well-known advantages that several species offer as experimental models. A remarkable feature of this genus that comprises more than two thousand species [21] is its diverse ecology: some species use fruits as breeding sites, others use flowers, mushrooms, sap fluxes, and fermenting cacti (reviewed in [22-25]). The adoption of decaying cacti as breeding sites occurred more than once in the evolutionary history of Drosophilidae [26, 27] and is considered a key innovation in the diversification and invasion of American deserts and arid lands by species of the Drosophila repleta group (hereafter the repleta group) [26]. Many species in this group are capable of developing in necrotic cactus tissues and feeding on cactophilic yeasts associated to the decaying process [28-35].

The repleta group comprises more than one hundred species [23, 36–39], however, only one of the more than twenty complete (or nearly complete) Drosophila mitogenomes assembled so far belongs to a species in this group (checked in GenBank, March 28, 2019), D. mojavensis (GenBank: BK006339.1). The latter, the first cactophilic fly to have a sequenced nuclear genome [40], is a member of the D. mulleri complex, an assemblage of species that belongs to the D. mulleri subgroup, one of the six species subgroups of the repleta group [37].

The D. buzzatii complex is the sister group of the D. mulleri complex [26]. It diversified in the Caribbean islands and South America, giving rise to the D. buzzatii cluster (hereafter the buzzatii cluster), and the D. martensis and D. stalkeri clusters [41]. The former is an ensemble of seven closely related species including D. antonietae [42], D. borborema [43], D. buzzatii [44], D. gouveai [42], D. koepferae [45], D. serido [43], and D. seriema [46]. All species are endemic to South America (Fig 1), except the semi-cosmopolitan D. buzzatii that reached a wide distribution following human mediated dispersion of prickly pears in the genus Opuntia (Caryophillales, Cactaceae) in historical times [35, 47, 48]. These species inhabit open areas of sub-Amazonian semidesert and desert regions of South America, where flies use necrotic cactus tissues as obligatory feeding and breeding resources [35, 49]. Regarding host plant use, D. buzzatii is an Opuntia specialist [31], considered an ancestral condition [26]. However, D. buzzatii has also been reared from necrotic columnar cacti [35]. The remaining species are mainly columnar dwellers although D. antonietae and D. serido can also use O. monacantha [49], while D. koepferae has also been recovered from decaying cladodes of O. monacantha, O. quimilo and O. sulphurea [31].

Fig 1

Geographical distribution of buzzatii cluster species modified from reference [50].

Species of the buzzatii cluster are almost indistinguishable in external morphology, however, differences in the morphology of the male intromittent organ (aedeagus) and polytene chromosome inversions provide clues to species identification (reviewed in [35, 48, 51]). The cluster has been divided into two groups based on aedeagus morphology, the first includes D. buzzatii and the remaining species compose the so-called Drosophila serido sibling set -serido sibling set from hereafter- [48]. In turn, analysis of polytene chromosomes revealed four informative paracentric inversions that define four main lineages: inversion 5g fixed in D. buzzatii, 2j in D. koepferae, 2x shared by D. antonietae and D. serido, and 2e shared by D. borborema, D. gouveai, and D. seriema [41, 52]. However, neither genital morphology nor chromosomal inversions are useful for inferring basal relationships within the cluster.

Pre-genomic phylogenetic studies based on a few molecular markers generated debate since different tree topologies were recovered depending on the molecular marker used. On one hand, the mitochondrial cytochrome oxidase I (COI) and the X-linked period gene supported the hypothesis of two main clades, one including D. buzzatii and D. koepferae and another comprising the remaining species [48, 53, 54]. On the other hand, trees based on a few nuclear and mitochondrial markers supported the hypothesis that D. koepferae was sister to the serido sibling set [26, 55]. To further complicate this issue, not all the same species were analyzed in these studies. In this vein, a recent genomic level study using a large transcriptomic dataset supports the placement of D. koepferae in the serido sibling set and D. buzzatii as sister to this set [50] similar to the results of [26]. However, phylogenetic relationships within the serido sibling set could not be ascertained despite the magnitude of the dataset employed by Hurtado and co-workers [50]. Thus, our aim is to shed light on the evolutionary relationships within the buzzatii cluster by means of a mitogenomic approach.

In this paper, we report the assembly of the complete mitogenomes of D. antonietae, D. borborema, D. buzzatii, D. seriema and two strains of D. koepferae, together with the corresponding gene annotations. Unfortunately, D. gouveai and D. serido, that inhabit Brazilian arid lands, could not be included because they are difficult to obtain and are not available in Drosophila repositories. We also present a mitogenomic analysis that defines a different picture of the relationships within the buzzatii cluster with respect to the results generated with nuclear genomic data. Finally, we discuss possible causes of the discordance between nuclear and mitochondrial datasets.

Material and methods

Species selection

The mitochondrial genomes of six isofemale lines of five species of the buzzatii cluster were assembled for the present study, for which NGS data were available. D. antonietae (Dato) was collected in March 2010 on Martín García Island (Buenos Aires Province, Argentina 34°10′42.12″S 58°15’23.15”W) by J. Hurtado and E. Hasson. D. borborema (Dbrb) was obtained from the Drosophila Species Stock Center (Stock Number: 15081–1281.01, University of California, San Diego, USA) and derived from collections performed in 1974 on Morro do Chapéu (Bahía State, Brazil 11°34′51.98″S 41°07′08.13″W) by M. Wasserman and R.H. Richardson. D. buzzatii (Dbuz) was collected in summer 2010 in Lavalle (Mendoza Province, Argentina 32°37′26.44″S 67°34′15.20″W) by J. Hurtado and E. Hasson. Two D. koepferae (Dkoe) strains; strain B from collections made in Bolivia in December 1982 by A. Fontdevila and A. Ruiz, and strain A collected in Vipos (Tucumán, Argentina 26°28′59″S 65°22′00″W) in February 2010 by J. Hurtado and E. Hasson. D. seriema (Dsei, strain D73C3) derived from collections on Cachoeira do Ferro Doido (Bahía State, Brazil 11°37' 40'' S 41°9' 2'' W) in June 1990 by G. Kuhn and F.M. Sene [56]. The stocks of D. antonietae, D. borborema, D. buzzatii and D. koepferae are available upon request. The rationale of including these D. koepferae strains is motivated by previous protein electrophoresis work showing higher genetic divergence between Bolivian and Argentinian populations than between conspecific populations in other species [45]. In addition, we also included four species of the subgenus Drosophila, for which assembled mitogenomes were available as outgroups in the phylogenetic analyses: D. grimshawi (GenBank: BK006341.1), D. littoralis (GenBank: NC_011596.1), D. virilis (GenBank: BK006340.1) and D. mojavensis (GenBank: BK006339.1).

In silico mtDNA reads extraction

Whole genome sequencing (WGS) and RNA-seq data for D. antonietae, D. borborema and both strains of D. koepferae were generated in our laboratory [34, 50]. For D. seriema and D. buzzatii, mitochondrial reads were retrieved from the Genome sequencing of D. seriema deposited in Sequence Read Archive database (SRA accession ID: ERX2037878) [56] and the D. buzzatii genome project (https://dbuz.uab.cat), respectively. For each species, mitochondrial reads were extracted from genomic and transcriptomic (when available) datasets. Bowtie2 version 2.2.6 [57] was first used with default parameters (end-to-end sensitive mode) to map reads to the mitochondrial genome of D. mojavensis, the closest relative of buzzatii cluster species available, as a reference. Next, only reads that correctly mapped to the reference genome were retained using Samtools version 1.8 [58]. Finally, mapped reads from genomic and transcriptomic datasets were combined to generate a set of only mitochondrial reads.

Mitochondrial reference genome assembly

It is well known that at least 25% of NGS reads are of mitochondrial origin [3]. Therefore, after the mapping process it is possible to attain a coverage ranging from 2000x to more than 20000x for mitogenomes. In order to avoid misassemblies caused by a large number of reads and given the difficulty of determining the coverage (and combination of reads) that recovers the complete mitochondrial genome, we split the reads into several datasets with different coverages by random sampling. Then, a three-step assembly procedure was adopted for these datasets based on recommendations of MITObim package version 1.8 [1]. In the first step, each dataset was employed to build a template by mapping its reads to the mitogenome of D. mojavensis using MIRA assembler [59]. In this way, several templates, based mostly on conserved regions, were built for each species. In the second step, entire mitogenomes were assembled by mapping the complete set of reads to the templates generated in the first step (coverage assembly), individually. This step was performed with the MITObim script and a maximum of ten mapping iterations. Finally, all the different coverage assemblies of the same species were aligned with Clustalw2 version 2.1 [60], and a consensus assembly was then generated considering a sequence representation threshold of 60% (nucleotide mostly represented at each position between assemblies), and not allowing gaps. De novo assemblies for each species, though more fragmented, were aligned to the assemblies obtained as described above and revealed the same gene order along the mitogenomes.

PCR amplification, sanger sequencing and consensus correction

Mitogenome assembly coverage averaged more than 20000x; however, three regions including parts of COI, NADH dehydrogenase subunit 6 (ND6) and large ribosomal RNA (rRNAL) genes showed low read representation in all species, producing miss-assemblies and fragmentation. These regions were PCR-amplified with GO taq Colorless Master Mix by Promega using primers designed for regions conserved across the six mitogenomes assembled in this study (data in S1 Text). PCR amplifications included an initial denaturation at 94°C for 90 s, followed by 25 cycles of denaturation at 94°C for 45 s, annealing at 62°C for 50 s, extension at 72°C for 1 min and a final 4 min extension. PCR fragments were sequenced in both directions on an ABI-3130xl (Genetic Analyzer). Sequences were analyzed and filtered using Mega X software [61] and, finally, merged with the assemblies.

Genome annotation and bioinformatic analyses

The six new assemblies were annotated with the MITOS web server (http://mitos.bioinf.uni-leipzig.de) [62] using the invertebrate mitochondrial genetic code and default parameter settings. The position and orientation of annotations were examined by mapping reads to mitogenomes with Bowtie2 [57] and visualization conducted with IGV ver. 2.4.10 [63]. In addition, nucleotide composition and codon usage were analyzed using MEGA X [61]. A homemade python package (available upon request) was developed to compute the number of pairwise nucleotide differences in the buzzatii cluster, and to visualize its variation along the mitogenomic alignment. Then we used the p-distance as a measured of nucleotide divergence, by dividing the number of nucleotide differences by the total number of nucleotides compared and by the number of pairwise comparisons [61]. Similar p-distance estimates were computed for the D. melanogaster subgroup with the aim to compare divergence patterns along the mitochondrial DNA in the buzzatii cluster with a deeply studied ensemble of species. To this end, one mitogenome of each one of the following species: D. melanogaster (KJ947872.2), D. erecta (BK006335.1), D. simulans (NC_005781), D. sechellia (NC_005780) and D. yakuba (NC_001322.1) were aligned, and nucleotide divergence estimates (p-distance) were obtained as described above. Synonymous (dS) and non-synonymous substitution rates (dN) were also estimated for each mitochondrial protein coding gene (PCG) using PAML 4.8 [64]. These estimates, as well as the ω ratio (dN/dS), were obtained separately for the buzzatii cluster and the melanogaster subgroup sequence alignments. Multiple sequence alignments of each coding gene were obtained with Clustalw2 ver. 2.1 [60].

Phylogenetic analyses

Phylogenetic analyses were conducted considering PCGs, ribosomal genes (rRNAs), transfer RNA genes (tRNAs) and intergenic regions (excluding the control region) of the 6 mitogenomes plus the sequences of the outgroups D. virilis, D. grimshawi, D. littoralis and D. mojavensis (see details in species selection section). An alignment of the ten mitogenomes was performed with Clustalw2 version 2.1 [60]. The flanking sequences that correspond to the control region and portions of the alignment showing abundant gaps were manually removed with Seaview ver. 4 [65]. The final alignment was used as input in PartitionFinder2 [66] to determine the best partition scheme and substitution models, considering separate loci and codon position (in PCGs), which were used in Bayesian Inference and Maximum Likelihood phylogenetic searches. In the Bayesian Inference approach executed with MrBayes ver. 3.2.2 [67], both substitution model and parameter estimates were unlinked. Then, two independent Markov Chain Monte Carlo (MCMC) were run for 30 million generations with three samplings every 1000 generations, for a total of 30,000 trees. Tracer ver. 1.7.1 [68] was used to assess the convergence of the chain mixing, where all parameters had effective sample sizes (ESS) > 200, and 25% of the trees were discarded as burn-in and the remaining trees were used to estimate a consensus tree and the posterior probability of each clade. The consensus tree was plotted and visualized with FigTree ver. 1.4.4 (https://github.com/rambaut/figtree/releases) [69]. Maximum Likelihood searches were performed in 2,000 independent runs using RAxML ver. 8.2.11 [70], applying the rapid hill climb algorithm and the GTR+GAMMA model, considering the partition scheme obtained with PartitionFinder2. Two thousand bootstrap replicates were run to obtain clade frequencies that were plotted onto the tree with highest likelihood. Tree and bootstrap values were visualized with FigTree ver. 1.4.4 [69]. Bayesian Inference searches for each PCG were individually performed to identify correlations with the topology recovered using the complete mitogenome. The GTR-GAMMA model together with the same parameters and evaluation detailed before were applied on each MCMC.

Divergence time estimation

Divergence times were estimated using the same methodology as in Hurtado et al. [50]. Four-fold degenerate third codon sites, i.e. putative neutral sites, of PCGs were extracted from the alignment and Bayesian Inference searches were analyzed with BEAST ver. 1.10.4 [71]. A strict clock was set using a prior for the mutation rate of 6.2x10-07 per year (standard deviation of 1.89x10-07), as was empirically estimated for mitochondrial DNA in Drosophila melanogaster [72]. In addition, a birth-death process with incomplete sampling and a time of 11.3 million years ago -MYA- (CI = ~9.34 to ~13) [26] to the root were defined as tree priors. Two MCMC were produced in 30 million generations with tree sampling every 1000 generations. Tracer [68] was used to evaluate chain convergence, discarding 10% of the all trees (burn-in). The information of the recovered trees was summarized in one tree applying LogCombiner and TreeAnnotator ver. 1.10.4 (available as part of the BEAST package), including the posterior probabilities of the branches, the age of the nodes, and the posterior estimates and HPD limits of the node heights. The target tree was visualized using FigTree [69]. Only D. mojavensis was included as an outgroup in this analysis to minimize problems of among-taxa rate variation given by the large divergence between the buzzatii cluster and the rest of the species already sequenced, together with the lack of time point calibrations and accurate mutation rates.

Results

Mitogenome characterization, nucleotide composition and codon usage

The length of the assembled mitogenomes varied from 14885 to 14899 bp among the six strains reported in this paper. Mitogenomes consisted of a conserved set of 37 genes, including 13 PCGs, 22 tRNAs and 2 rRNAs genes, with order and orientation identical to D. mojavensis. Several short non-coding intergenic regions were also found. Twenty-three genes were found on the heavy strand (+) and fourteen on the light strand (-). Detailed statistics about metrics and composition of the mitogenomes are shown in Table 1.

Table 1

Composition of mitochondrial elements in the species assemblies of the Drosophila buzzatii cluster.

	Species assembly
Composition	D. antonietae	D. borborema	D. buzzatii	D. koepferaeB	D. koepferaeA	D. seriema	D. mojavensis
Total length	14885	14889	14889	14892	14891	14891	14904
GC (%)	23.36	23.22	23.60	23.20	23.20	23.26	23.54
N's (%)	0.00	0.00	0.00	0.00	0.00	0.00	0.00
Intergenic (%)	2.76	2.71	2.64	2.80	2.77	2.90	0.93
tRNAs (%)	9.82	9.81	9.80	9.81	9.81	9.80	9.78
PCGs (%)	72.98	73.04	73.15	72.95	72.97	72.85	75.15
rRNAs (%)	14.44	14.45	14.41	14.44	14.44	14.44	14.14

GC, Guanine+Cytosine; N's: undetermined nucleotides; Intergenic, regions between genes (non-coding); tRNAs, transfer RNA; PCGs, Protein Coding Genes; rRNAs, ribosomal RNA.

Overall nucleotide composition in PCGs ranged between 37.6–37.8% A, 37.2–37.9% T, 10.2–10.4% G, and 14.1–14.7% C. The thirteen PCGs were AT-biased as in the entire mitogenome, and the codon usage bias in each gene was greater than 0.50. The most frequently used codons were UUA (Leu), AUU (Ile), and UUU (Phe) in all cases. Codon usage information for each species is shown in Table in S1 Table.

Genetic diversity among mitogenomes

Patterns of divergence (p-distance) along the entire mitogenomes were, overall, very similar for both the buzzatii cluster and the melanogaster subgroup (Fig 2). However, overall divergence was higher in the melanogaster subgroup than in the buzzatii cluster. This difference between both ensembles is probably due to depth of divergence, the melanogaster subgroup comprises pairs of highly diverged species such as D. melanogaster and D. yakuba that split more than 10 MYA [73]. Despite the general similarity in the patterns of divergence, a substantial difference was found in the region encompassing COIII, tRNA-G and ND3 genes. In this region, from position 5000 to 6000, nucleotide diversity was the highest in the melanogaster subgroup, showing an apparent increase represented by two high peaks absent in the buzzatii cluster. Considering genetic divergence within the buzzatii cluster (Fig 3), the lowest value of average pairwise nucleotide divergence was observed for the pair D. borborema and D. seriema (p = 1.91x10-03), while between D. seriema and D. buzzatii divergence was an order of magnitude larger (p = 2.73x10-02). Divergence between D. koepferae strains was surprisingly high (7.14x10-03). The complete set of divergence estimates in the buzzatii cluster is reported in Table in S2 Table. Substitution rates in synonymous (dS) and non-synonymous (dN) sites, and the ω ratios varied among PCGs (Table 2). The dN/dS ratio (ω) varied from 0.003 to 0.060 among PCGs in the buzzatii cluster. The range in the melanogaster subgroup was similar, but with a lower upper bound (0.003–0.018). Two genes, ATP8 and ND2, appear as outliers in the buzzatii cluster (ATP8 and ND2), which apart from these two loci, had lower divergence values than in the melanogaster subgroup. The elevated dN/dS values observed for ATP8 and ND2 in the buzzatii cluster (Table 2) are probably real since the true value of dS may be slightly underestimated due to multiple hits in the melanogaster subgroup (see above), leading to a slight overestimation of the dN/dS ratio. In any case, these results suggest that purifying selection imposes strong constraints in the evolution of mitochondrial genes.

Fig 2

Nucleotide diversity variation along the mitogenome, estimated for a sliding window of 500bp with an overlap of 100bp.

p-distance values for species belonging the buzzatii cluster and the melanogaster subgroup are represented independently.

Fig 3

Pairwise comparison of nucleotide diversity between species belonging the buzzatii cluster.

The official FlyBase abbreviations for Drosophila species names are shown.

Table 2

Estimates of non-synonymous (dN) and synonymous (dS) substitutions and their ratio (ω) among species of the buzzatii cluster and the melanogaster subgroup.

	buzzatii cluster			melanogaster subgroup
PCG	ω	dN	dS	ω	dN	dS
ATP6	0.005	0.007	1.390	0.014	0.034	2.364
ATP8	0.060	0.031	0.506	0.014	0.071	4.892
CytB	0.005	0.008	1.463	0.010	0.027	2.556
COI	0.003	0.003	1.211	0.003	0.006	2.383
COII	0.005	0.005	1.058	0.008	0.018	2.295
COIII	0.006	0.008	1.277	0.012	0.024	1.975
ND1	0.003	0.012	4.332	0.006	0.028	4.490
ND2	0.036	0.040	1.102	0.016	0.039	2.368
ND3	0.009	0.034	3.657	0.016	0.051	3.085
ND4l	0.008	0.009	1.158	0.001	0.013	8.026
ND4	0.011	0.013	1.253	0.009	0.040	4.306
ND5	0.007	0.018	2.564	0.013	0.063	4.764
ND6	0.012	0.027	2.231	0.018	0.082	4.534

The sequences of the 13 PCGs, 22 tRNAs genes, 2 rRNAs genes, and intergenic regions were included in the alignment. Total length of the final matrix encompassing the ten mitogenomes was 15044 characters, from which 1950 were informative sites, 11583 conserved, and 1422 were singletons. Both Maximum Likelihood and Bayesian Inference phylogenetic analyses recovered the same highly supported topology that confirmed the monophyly of the buzzatii cluster (Fig 4). Two main clades can be observed in the tree, one including both D. koepferae strains as sisters to D. buzzatii, and the second, comprising D. antonietae as sister species of the sub-clade formed by D. borborema and D. seriema. The species selected as outgroups were placed as expected, with D. mojavensis as the closest relative of the buzzatii cluster. We also performed a gene tree analysis using all PCGs (S1 Fig). We could only obtain trees for seven genes out of the thirteen PCGs, given the lack of informative sites in the alignments of ATP8, ATP6, ND3, ND4l, COII and COIII. Only two (ND1 and ND5) out of the seven recovered gene trees showed the same topology as the complete mitogenome, while the remaining genes produced three (different) topologies. Trees obtained with CytB and ND4 showed D. buzzatii as sister to the serido sibling set which included D. koepferae. COI and ND2 retrieved trees where D. buzzatii and D. koepferae exchanged positions in the tree, placing D. koepferae as the species closest to the putative ancestor of the cluster. ND6 recovered two clades where D. antonietae was the sister of D. buzzatii and D. koepferae (both strains) in one clade, and the pair D. borborema-D. seriema composed the other.

Fig 4

Phylogenetics hypotheses for the buzzatii cluster based on the entire sequence of the mitogenome (control region not included).

Tree topology recovered by both Maximum Likelihood and Bayesian Inference searches. Bootstrap values and posterior probabilities are indicated at each node.

Divergence times

PCGs contained 1201 4-fold degenerate sites in the mitogenomes of the buzzatii cluster strains assembled in this study. The tree obtained in the divergence time estimation analysis (Fig 5) was topologically identical to the trees obtained in the phylogenetic analyses using complete mitogenomes (see Fig 4). Divergence time estimations showed that the buzzatii cluster diverged in the Early Pleistocene, 2.11 MYA, and the split with the D. mojavensis common ancestor occurred 10.63 MYA in the Miocene. Our results also indicated that the clade containing D. antonietae, D. borborema and D. seriema, is younger than the clade composed by D. buzzatii and D. koepferae. In addition, the split between D. borborema and D. seriema is quite recent, about ~50,000 years ago, in the Late Pleistocene, even more recent than the split of D. koepferae strains that diverged ~310,000 years ago, in the Middle Pleistocene.

Fig 5

Divergence times for the buzzatii cluster drawn on a Bayesian inference tree.

Numbers on each node are the time estimates. Blue bars represent the 95% confidence intervals of estimates.

Discussion

The six newly assembled mitochondrial genomes of five cactophilic species of the buzzatii cluster share molecular features with animal mitochondrial genomes sequenced so far [74]. All assembled mitogenomes contain the same set of genes usually found in animal mitochondrial genomes. Gene order and orientation, as well as the distribution of genes on the heavy and light strands were identical to the mitogenome of the closest relative D. mojavensis and other drosophilids [9]. Analysis of overall nucleotide composition of mitogenomes and PCGs revealed the typical AT-bias found in Drosophila mitogenomes. Codon usage is highly biased suggesting that synonymous sites cannot be considered strictly neutral and that some sort of natural selection for translational accuracy governs codon usage [75].

Phylogenetic analyses based either on complete mitogenomes or four-fold degenerate sites (for divergence time estimations), retrieved a well supported tree, suggesting that the cluster is composed by two main clades, one including D. buzzatii and D. koepferae (both strains) and another comprising D. antonietae, D. seriema and D. borborema. Though our present results are consistent with previous work based on single mitochondrial genes [53, 76], they should be considered with caution since we only included a single inbred line as representative of each species, except for D. koepferae. Moreover, phylogenetic relationships depicted by mitogenomes do not agree with phylogenetic studies based on both a small set of nuclear and mitochondrial genes [26] and a large set of nuclear genes -see below- [50]. Interestingly, the topology showing these two clades was only recovered in two (ND1 and ND5) out of seven trees based on individual PCGs, and the remaining gene trees produced either a novel topology or a topology consistent with the phylogeny reported in Hurtado et al. [50].

The lack of recombination causes mitochondrial DNA to be inherited as a unit, so trees recovered with individual mitochondrial genes are expected to share the same topology and to be consistent with the trees obtained with complete mitogenomes. On the other hand, our results suggest that individual genes not only produce different topologies but also a poor resolution of phylogenetic relationships. Such inconsistencies between complete mitogenomes and gene trees in phylogenetic estimation may result from inaccurate reconstruction or from real differences among gene trees. The first possible explanation is the simple fact that numbers of informative sites within a single locus are insufficient to accurately estimate phylogenetic relationships, particularly in groups of recently diverged species: overall support for gene trees was poorer than for the tree based on complete mitogenomes. Second, heterogeneity in evolutionary rates among genes and/or differences in selective constraints along the mitogenome can also account for these inconsistencies [17, 77–79] consistent with our observations of substantial variation in synonymous and nonsynonymous rates as well as ω ratios across PCGs. In addition, variation among oxidative phosphorylation complexes in the buzzatii cluster was high. The ND complex was, on average, less constrained than the ATP complex, cytochrome b (CytB) and cytochrome oxidase complex (COI, COII & COIII), consistent with results reported in the melanogaster subgroup [9, 80]. Another factor that may lead to biased tree construction, particularly relevant for mitochondrial genes characterized by high substitution rates, is substitutional saturation [81]. A priori, saturation should not be problematic in recently diverged species, like the buzzatii cluster, however, saturation may be problematic in the estimation of divergence relative to the outgroup and, thus, for phylogenetic inference. The closest outgroup to the buzzatii cluster employed in our study was the mulleri complex species D. mojavensis. Available evidence suggest that these complexes diverged ~10 MYA [26], but possibly more recently, ca 5.5 MYA [50] suggesting that substitution saturation may lead to inaccurate phylogenetic reconstruction.

In this context, a recent report investigating the effect of using individual genes, subsets of genes, complete mitogenomes and different partitioning schemes on tree topology suggested a framework to interpret the results of mitogenomic phylogenetic studies [11]. The authors concluded that trees obtained with complete mitogenomes reach the highest phylogenetic performance and reliability than single genes or subsets of genes. Therefore, we consider that phylogenetic relationships inferred from complete mitogenomes reflect the evolutionary history of, at least, mitogenomes.

The phylogenetic relationships depicted by our mitogenomic approach are incongruent with a recent study based on transcriptomic data [50]. Based on a concatenated matrix of 813 kb uncovering 761 gene regions, the authors obtained a well-supported topology in which D. koepferae appears phylogenetically closer to D. antonietae and D. borborema than to D. buzzatii, placing D. buzzatii alone as sister to the rest of the cluster. This topology agrees with male genital morphology, cytological and molecular phylogenetic evidence [26, 35, 55]. Nevertheless, the pattern of cladogenesis of the trio D. koepferae-D. borborema-D. antonietae could not be fully elucidated since a nuclear gene tree analysis yielded ambiguous results. As a matter of fact, the analysis of the 761 gene trees reported showed that about one third of the genes supported each one of the three possible topologies for the trio D. koepferae-D. antonietae-D. borborema indicating a hard polytomy [50]. In contrast, the early separation of D. buzzatii from the serido sibling set is supported by 97% of the gene trees and, surprisingly, none of the gene trees recovered the clade including D. buzzatii and D. koepferae as the sister group of the clade involving D. antonietae and D. borborema [50] as suggested by the present mitogenomic approach. Such mitonuclear discordance has been reported in several animal species. A recent review lists several examples in animals [82]. Likewise, the literature in this respect is abundant in the genus Drosophila. Well-known cases are D. pseudoobscura and D. persimilis [83]; D. santomea and D. yakuba [84]; and D. simulans and D. mauritiana [85]. Mitonuclear discordance may be caused by incomplete lineage sorting (ILS) and/or introgressive hybridization. These two factors do not equally affect mitochondrial and nuclear genomes because ILS is more likely for nuclear genes, especially when the ancestral effective population size of recently diverged species was large [86, 87]. Introgressive hybridization is expected to be prevalent in mitochondrial genomes given its lower effective population size [88]. If we accept that the topology based on nuclear genes is representative of the species-history (see also [49]), the greater similarity between D. buzzatii and D. koepferae mitogenomes is suggestive of gene flow between these largely sympatric species [35]. Thus, we suggest that D. buzzatii and D. koepferae lineages initially separated but then exchanged genes via fertile F1 females (males were likely sterile as expected due to Haldane’s Rule) before finally separating less than 1.5 MYA. Not only the more recent mitogenomic ancestry is suggestive of gene exchange, also traces of introgressive hybridization can still be detected in nuclear genomes [50].

In fact, phylogenetic, population genetic and experimental hybridization studies suggest a significant role of introgression in the evolutionary history of the buzzatii cluster. Phylogeographic studies revealed discordances between mitochondrial markers and genital morphology in areas of sympatry between species [53]. Likewise, interspecific gene flow has been invoked to account for shared nucleotide polymorphisms in nuclear genes in D. buzzatii and D. koepferae that cannot be accounted by ILS [89, 90]. Moreover, experimental hybridization studies have shown that several species of the buzzatii cluster can be successfully crossed, producing fertile hybrid females that can be backcrossed to both parental species. Interestingly, D. koepferae can be successfully crossed with D. antonietae, D. borborema, D. buzzatii and D. serido [45, 48, 91–95].

Our estimates of divergence times are in conflict with most previous studies. In general, previous estimates, based on individual genes or a few genes (either mitochondrial or nuclear) suggested an older origin of the cluster and deeper splitting times within the cluster when compared to the estimates based on transcriptomes and mitogenomes. In effect, Gómez & Hasson [89] and Oliveira et al. [26] dated the split of D. buzzatii from the remaining species of the cluster at 4.6 MYA, respectively, whereas Manfrin et al.’s estimates [48] were even older, from 3 to 12 MYA for the most recent to the more ancient split. In contrast, putting apart the divergence time of the clade D. buzzatii-D. koepferae for the reasons discussed above, the radiation of the remaining three species seems to be extremely recent, less than 1 MYA using mitochondrial genomes (Fig 5), which are similar to estimates based on transcriptomes [50]. However, it is worth mentioning that divergence times estimated in the present paper and by Hurtado et al. [50] may be biased downwards since both are based on empirical mutation rates for nuclear and mitochondrial genes, respectively, calculated for over 200 generations in D. melanogaster [72]. Thus, these results should be interpreted with caution in the light of evidence suggesting not only the time-dependence of molecular evolutionary rates, but also that mutation rates obtained using pedigrees and laboratory mutation-accumulation lines often exceed long-term substitution rates by an order of magnitude or more [79].

Even though divergence times estimates obtained in this study cannot be entirely compared to assessments based on nuclear genomic data and individual nuclear genes, given the uncertainty of tree topology, they concur in that species of the buzzatii cluster apparently emerged during the Late Pleistocene in association with Quaternary climate fluctuations [49, 50, 76]. Moreover, in view of the obligate ecological association between buzzatii cluster species and cacti, the so-called Pleistocene “refuge hypothesis” is a suitable explanation for the diversification in this group in active cladogenesis. This hypothesis suggests that Pleistocene glacial cycles successively generated isolated patches of similar habitats across which populations may have diverged into species [96, 97].

Available paleo-climatic evidence, consistent with the Pleistocene “refuge hypothesis”, can also account for the relatively deep intraspecific divergence between Bolivian and Argentinian D. koepferae strains. Because Quaternary topographical patterns in the Central Andes have remained unchanged for the last 2–3 MYA, a plausible explanation for this late Pleistocene vicariant event is related with glacial-interglacial cycles [98]. Although the validity of the Pleistocene “refuge hypothesis” is controversial (cf. [99]) and few studies addressed specific hypotheses on how the Quaternary glacial-interglacial cycles impacted species diversification [100], our divergence time estimates between Bolivian and Argentinian D. koepferae suggest a role of climatic oscillations as a factor of ecogeographical isolation in the Central Andes during the Pleistocene. Moreover, paleo-climatological evidence suggest that the area inhabited by D. koepferae has been exposed to substantial climatic variations on timescales of 103–105 years related with glacial-interglacial cycles. Thus, Andean north-south exchanges may have been alternately favored or disfavored by these Quaternary climatic oscillations. In fact, the estimated age of the vicariant event between the D. koepferae strains is tantalizingly coincident with the coldest phase of the Marine Isotopic Stage (MIS) 10, which corresponds to a glacial period that ended about 337,000 years ago [101]. The coldest period of the MIS 10, recorded in global air and sea surface temperature and also the lowest atmospheric CO2 levels, occurred ca 355,000 years, well within the confidence interval of our divergence time estimated between D. koepferae strains. On a global scale, glacial periods are primarily reflected in a lowering of air temperature but also in altered patterns of precipitation in the both sides of the Central Andes [102] which were in turn the main drivers of vegetation changes [103] including the appearance of South American columnar cacti [104]. Besides the impact on air temperature, periods of ice advance in the Central Andes generally were periods of negative water balance in the Pacific coastal regions west to the Central Andes [105], and a positive water balance in the Central Andes, as evidenced by deeper and fresher conditions in Lake Titicaca [106] (see S2 Fig). Thus, during the colder and wetter phases of the MIS 10 in the Central Andes, species distributions may have suffered a general contraction towards the southern and northern lowland warmer refugia between 1000–2000 m, whereas a general worsening condition occurred in higher western elevations. Northern and southern refugia were probably separated by a gap of low suitability represented by the steep gradient of the eastern flank of Eastern Andes between 22–24°S, which represents today a region of strong W-E precipitation gradient. The MIS 10 glacial cycle has a particular structure since it does not have a pronounced interstadial (relative warmer) conditions in the mid-cycle [107], providing a prolonged, effective “soft” dispersal barrier that affected the distribution of D. koepferae.

Finally, future analyses including the mitogenomes of the other Brazilian species D. gouveai and D. serido and several mitogenomes of each species are needed to achieve a more complete understanding of the evolutionary history of the cluster. Such comparative analysis including the complete mitogenomes of all buzzatii cluster species will help to disentangle the intricate relationships in this group.

Set of primers used to sequence each gene.

(DOCX)

Click here for additional data file.

Codon usage for each mitogenome of the buzzatii cluster species.

(XLSX)

Click here for additional data file.

Genetic divergence among species of the buzzatii cluster.

Estimates are shown for each pairwise comparison between species.

(XLSX)

Click here for additional data file.

Phylogenetic hypotheses for the buzzatii cluster species recovered by each mitochondrial gene using Bayesian Inference searches.

(TIF)

Click here for additional data file.

Paleoclimatic records of the last 500,000 years.

Ages in the top are indicated as 103 years (kyrs). Gradated shading area indicates divergence age estimates. Marine Isotope Stages (MIS) are labeled according to Lisiecki and Raymo [101]. Shaded vertical areas correspond to glacial periods whereas white areas correspond to interglacials or interstadials. Glacial periods correspond to cold and dry conditions in the western slopes of the Western Andes, and cold and wetter conditions in the eastern slopes of the Eastern Andes and the Altiplano. A. Globally-averaged surface air temperature anomaly reconstructed from proxy and model data for the last eight glacial cycles [108]. B. CO2 concentration based on Vostok Ice Core data [109]. C. Iron accumulation rates (AR Fe) reflecting changes in terrigenous sediment input to ODP Site 1239D, Equatorial Pacific [110]. D. % of CaCO3from Site LT01-2B indicating changes in water balance at Lake Titicaca Basin, Bolivia (modified from [106]).

(TIF)

Click here for additional data file.

30 Aug 2019

PONE-D-19-20340

What does mitogenomics tell us about the evolutionary history of the Drosophila buzzatii cluster (repleta group)

PLOS ONE

Dear Mr. Moreyra,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Two external referees and me have now evaluated your paper. Both reviewers raised some very important issues that need to be addressed in a revision. I re-read the MS and found a large number of grammar and translation issues that must be corrected satisfactorialy in a revised MS. I have uploaded a marked up manuscript to aid you in the revision process. I have made many comments that need to be addressed in this marked up manuscript, but I may have missed some that need fixing. Please proofread carefully.

A few other points need explanation, such as why some in-group species were not included in this analysis, and some cautionary explanation about the use of single inbred lines to represent these species, except for D. koepferae.  I'd guess that more geographically disparate samples would complicate this analysis and change your conclusions as well. Further, using a "total evidence" approach, why not use all available genetic data including inversions to build these trees?  I don't think its surprising that mtDNA yields different results than nuclear markers. However, the acquisition of more mitogenomes for these species is a valuable exercise.

We would appreciate receiving your revised manuscript by Oct 14 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

William J. Etges

Academic Editor

PLOS ONE

Journal Requirements:

1. When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. In your Methods section, please provide additional location information of the collection sites, including geographic coordinates for the data set if available.

3. In your Methods section, please provide additional information regarding the permits you obtained for the work. Please ensure you have included the full name of the authority that approved the field site access and, if no permits were required, a brief statement explaining why.

4. We note that Figure 1 in your submission contain map/satellite images which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.

We require you to either (1) present written permission from the copyright holder to publish these figures specifically under the CC BY 4.0 license, or (2) remove the figures from your submission:

You may seek permission from the original copyright holder of Figure 1 to publish the content specifically under the CC BY 4.0 license.

We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:

“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”

Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.

In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”

If you are unable to obtain permission from the original copyright holder to publish these figures under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.The following resources for replacing copyrighted map figures may be helpful:

USGS National Map Viewer (public domain): http://viewer.nationalmap.gov/viewer/

The Gateway to Astronaut Photography of Earth (public domain): http://eol.jsc.nasa.gov/sseop/clickmap/

Maps at the CIA (public domain): https://www.cia.gov/library/publications/the-world-factbook/index.html and https://www.cia.gov/library/publications/cia-maps-publications/index.html

NASA Earth Observatory (public domain): http://earthobservatory.nasa.gov/

Landsat: http://landsat.visibleearth.nasa.gov/

USGS EROS (Earth Resources Observatory and Science (EROS) Center) (public domain): http://eros.usgs.gov/#

Natural Earth (public domain): http://www.naturalearthdata.com/

5. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

6. We note that you have included the phrase “data not shown” in your manuscript. Unfortunately, this does not meet our data sharing requirements. PLOS does not permit references to inaccessible data. We require that authors provide all relevant data within the paper, Supporting Information files, or in an acceptable, public repository. Please add a citation to support this phrase or upload the data that corresponds with these findings to a stable repository (such as Figshare or Dryad) and provide and URLs, DOIs, or accession numbers that may be used to access these data. Or, if the data are not a core part of the research being presented in your study, we ask that you remove the phrase that refers to these data.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: Comments to MS by Moreyra et al. PLOS ONE

General comments

This paper reports the sequencing, assembly and annotation of six mitochondrial genomes of five Drosophila repleta group species (subgenus Drosophila). This is useful information because none of these mitogenomes were previously available. The authors analyze these genomes to characterize patterns of divergence along the mitogenome, to elucidate the phylogenetic relationships between the five species and to estimate divergence times among them. It seems to me a straightforward, well written and useful paper. My comments are relatively minor.

Specific suggestions

Please, when necessary (e.g. Figure 3) use the official FlyBase abbreviations for Drosophila species names.

Dato = D. antonietae

Dbrb = D. borborema

Dbuz = D. buzzatii

Dkoe = D. koepferae

Dsei = D. seriema

To refer to the two strains of D. koepferae you should use always the same name: either A and B, or 11 and 7.1. Currently there is a mixture, compare Figure 3 and Figure 4.

Page 2, lines 29-30. I suggest to use the alphabetical order for species names (as in Table 1): D. antonietae, D. borborema, D, buzzatii, D. koepferae (two stocks) and D. seriema.

Figure 3. I don’t see the logical in the order of the pairwise comparisons. Perhaps you could use the alphabetical order (in agreement with previous comment):

Dato – Dbrb

Dato – Dbuz

Dato – Dkoe A

Dato – Dkoe B

Dato – Dsei

Dbrb – Dbuz

Dbrb – Dkoe A

Dbrb – Dkoe B

Dbrb – Dsei

Dbuz – Dkoe A

Dbuz – Dkoe B

Dbuz – Dsei

Dkoe A – Dkoe B

Dkoe A – Dsei

Dkoe B – Dsei

Pages 6-7, lines 144-155. It is not clear to me whether the authors produced a “de novo” assembly of the six mitogenomes or they used D. mojavensis as a template. This is important in order to conclude (page 10) that PCGs in the six mitogenomes show the same order and orientation than in D. mojavensis. Please clarify this issue.

Page 11, lines 246-247. Please give more details on how pairwise nucleotide diversity (�) was calculated. Because you have six mitogenomes in the buzzatii cluster, there is a total of 15 pairwise comparisons. Is the value reported in Figure 2 the average of the 15 pairwise comparisons? If this is the case, this value includes both interspecific and intraspecific (two stocks of D. koepferae) comparisons. Clarify.

Minor corrections

Page 2, line 33. The monophyly of that the D. buzzatii cluster. Delete “that”

Page 5, line 108. Insert “cluster” after buzzatii

Page 6, line 134. Split the sentence:

… unpublished data). For D. seriema and D. buzzatii, mitochondrial reads were retrieved …

Page 13, line 283. We could only obtain trees for ….

Page 13, line 287. Allocated -> placed

Page 13, line 294. Legend of Fig 4. Bootstrap values and posterior probabilities are indicated at each node.

Page 15, line 334. Thus -> on the other hand,

Heading of Table S1. Set of primers used to sequence each gene.

Reviewer #2: Moreyra et al. report on the assembly and phylogenetic analysis of mtDNAs from members of the Drosophila buzzatii group. Only one of the members of this clade had prior mtDNA assembly data, so this is a contribution to the genetics of this group of species. Using available genomic sequence data, the mtDNAs were assembled using bioinformatics tools, and gaps in coverage were validated using direct Sanger sequencing of amplified regions. The paper presents some straightforward molecular evolution and phylogenetic analyses that are consistent with other Drosophila mitogenomic analyses, but reveal some differences with the melanogaster group. The data are concordant with various biogeographic scenarios, and reveal some discordance with nuclear phylogenetic analyses. Below are a few comments that may improve the manuscript.

Line 33. Delete extra ‘that’

Line 246: A terminology point: ‘pi’ for pairwise diversity usually refers to polymorphism within populations, and estimates this value from a random sample of alleles. The use of pi from a sample of species has a somewhat different interpretation given that taxon sampling is different from a random sample of alleles from a breeding population. The use of a pairwise metric is fine, but some clarification should be made in these different contexts.

Line 248: The pi values are larger in the melanogaster subgroup than in the buzzatii cluster. Related to the comment above, if these were equivalent random samples from two different idealized populations the difference would imply different population genetic forces, but since the differences in divergence times and distribution of taxa in the two subgroups / clusters being compared are not the same, the pi metric may simply reflect those historical contingencies. The metric pi is sensitive to the frequencies of alleles in a population, so more closely related mtDNAs will reduce pi, or more deeper-diverging mtDNAs will inflate it. So, some comment should be made to this effect to provide context for the larger pi values in the melanogaster subgroup. Notably, that group has three very diverged mtDNAs within the D. simulans clade that may inflate pi, which is not seen in the buzzatii sample. Clarify that the Dmel vs. Dbuzz difference is probably due to depth of divergence, not something distinct about rates of evolution. The use of a pairwise metric is fine, but its interpretation is complicated by the difference between a taxon sample vs. an allele sample from a breeding population.

Lines 261-263: The higher dN/dS values for two genes in Dbuzz vs. Dmel data (Table2), may imply a shift in function. It could be noted that since the dS values (and pi as per points above) is higher in Dmel group, the true value of dS may be slightly underestimated due to multiple hits, leading to a slight overestimation of the dN/dS ratio in Dmel. This means that the elevated values of dN/dS for the ND and ATP genes in the Dbuzz group relative to the same Dmel genes are probably real.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: David Rand

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: Moreyra_et_al_2019_WJE.docx

Click here for additional data file.

18 Sep 2019

Editorial Board

Dr. William J. Edges

Academic Editorial

Dear Dr. Edges,

First, we would like to thank the subject editor Dr. W.J. Etges and the two reviewers, an anonymous one and Dr. D. Rand, for helpful and insightful comments that have definitely contributed to improve manuscript entitled “What does mitogenomics tell us about the evolutionary history of the Drosophila buzzatii cluster (repleta group)?” by Nicolás Moreyra et al.

We revised the manuscript addressing all the issues raised by the academic editor and both reviewers. Below, we transcribe each one of the comments, followed by our response (Author´s response), which outlines the changes made to the manuscript in order to reflect the suggestions made in each case.

Responses to the Academic Editor

Two external referees and me have now evaluated your paper. Both reviewers raised some very important issues that need to be addressed in a revision. I re-read the MS and found a large number of grammar and translation issues that must be corrected satisfactorialy in a revised MS. I have uploaded a marked up manuscript to aid you in the revision process. I have made many comments that need to be addressed in this marked up manuscript, but I may have missed some that need fixing. Please proofread carefully.

Author’s response: all comments and suggestions in the manuscript marked up by the editor were taken into account and incorporated in the revised version.

A few other points need explanation, such as why some in-group species were not included in this analysis, and some cautionary explanation about the use of single inbred lines to represent these species, except for D. koepferae. I'd guess that more geographically disparate samples would complicate this analysis and change your conclusions as well.

Author’s response: D. gouveai and D. serido are species that inhabit Brazilian arid lands, and could not be included in this study because of the unavailability of stocks in Drosophila repositories for these species. This is detailed in the revised manuscript. We included a sentence in the discussion section providing a cautionary note related to the use of single lines as representative of species, except for Dkoe, which was represented by two highly differentiated lines.

Further, using a "total evidence" approach, why not use all available genetic data including inversions to build these trees? I don't think its surprising that mtDNA yields different results than nuclear markers. However, the acquisition of more mitogenomes for these species is a valuable exercise.

Author’s response: We did not included other genetic data to build the trees because in a recent phylogenetic study based on a large nuclear genes dataset, considered other sources of information, like morphology and inversions. The aim of the present study was to compare the results of our mitogenomic approach with the study based on transcriptomic data.

Responses to reviewers

Reviewer #1:

Specific suggestions.

Please, when necessary (e.g. Figure 3) use the official FlyBase abbreviations for Drosophila species names.

Dato = D. antonietae

Dbrb = D. borborema

Dbuz = D. buzzatii

Dkoe = D. koepferae

Dsei = D. seriema

Author’s response: we adopted the official FlyBase abbreviations for the species names.

To refer to the two strains of D. koepferae you should use always the same name: either A and B, or 11 and 7.1. Currently there is a mixture, compare Figure 3 and Figure 4.

Author’s response: we corrected references to strains of D. koepferae to A and B (instead of 11 and 7.1) in the text, Figure 3, Figure 4 and Figure 5 (as well as in S1 Fig).

Page 2, lines 29-30. I suggest to use the alphabetical order for species names (as in Table 1): D. antonietae, D. borborema, D, buzzatii, D. koepferae (two stocks) and D. seriema.

Figure 3. I don’t see the logical in the order of the pairwise comparisons. Perhaps you could use the alphabetical order (in agreement with previous comment):

Dato – Dbrb

Dato – Dbuz

Dato – Dkoe A

Dato – Dkoe B

Dato – Dsei

Dbrb – Dbuz

Dbrb – Dkoe A

Dbrb – Dkoe B

Dbrb – Dsei

Dbuz – Dkoe A

Dbuz – Dkoe B

Dbuz – Dsei

Dkoe A – Dkoe B

Dkoe A – Dsei

Dkoe B – Dsei

Author’s response: We change the order of the species names in Table 1 to alphabetical order. However, we included D. mojavensis in the last column (avoiding the sorted list) as a comparative external reference.

Pages 6-7, lines 144-155. It is not clear to me whether the authors produced a “de novo” assembly of the six mitogenomes or they used D. mojavensis as a template. This is important in order to conclude (page 10) that PCGs in the six mitogenomes show the same order and orientation than in D. mojavensis. Please clarify this issue.

Author’s response: We produced a reference assembly only for the conserved regions between D. mojavensis and each species belonging the buzzatii cluster, creating several templates (based on different subsets of reads) for each species. After that, the complete set of mitochondrial reads (of each species) was mapped to the corresponding templates to generate several genome assemblies that were aligned to create a consensus per species. We also performed de novo assemblies for each species (line 166) that were more fragmented, but happily they revealed the same gene order along the mitogenomes. These details are clarified in the revised manuscript.

Page 11, lines 246-247. Please give more details on how pairwise nucleotide diversity (pi) was calculated. Because you have six mitogenomes in the buzzatii cluster, there is a total of 15 pairwise comparisons. Is the value reported in Figure 2 the average of the 15 pairwise comparisons? If this is the case, this value includes both interspecific and intraspecific (two stocks of D. koepferae) comparisons. Clarify.

Author’s response: we estimated nucleotide divergence by calculating the p-distance (analogous of pi), which is obtained by dividing the number of nucleotide differences between species pairs by the total number of nucleotides compared and by the number of pairs of sequences being compared (six in our study). We decided to independently include stocks of D. koepferae given the significant divergence between them, which is even greater than between D. borborema and D. seriema. In order to avoid misunderstandings regarding the parameter pi we decided to refers the measure of divergence as p-distance (instead of pi since this parameter usually refers to polymorphisms within populations).

Minor corrections

Page 2, line 33. The monophyly of that the D. buzzatii cluster. Delete “that”

Page 5, line 108. Insert “cluster” after buzzatii

Page 6, line 134. Split the sentence:

… unpublished data). For D. seriema and D. buzzatii, mitochondrial reads were retrieved …

Page 13, line 283. We could only obtain trees for ….

Page 13, line 287. Allocated -> placed

Page 13, line 294. Legend of Fig 4. Bootstrap values and posterior probabilities are indicated at each node.

Page 15, line 334. Thus -> on the other hand,

Heading of Table S1. Set of primers used to sequence each gene.

Author’s response: all corrections made

Reviewer #2:

Line 33. Delete extra ‘that’

Author’s response: corrected

Line 246: A terminology point: ‘pi’ for pairwise diversity usually refers to polymorphism within populations, and estimates this value from a random sample of alleles. The use of pi from a sample of species has a somewhat different interpretation given that taxon sampling is different from a random sample of alleles from a breeding population. The use of a pairwise metric is fine, but some clarification should be made in these different contexts.

Author’s response: see above in responses to reviewer 1.

Line 248: The pi values are larger in the melanogaster subgroup than in the buzzatii cluster. Related to the comment above, if these were equivalent random samples from two different idealized populations the difference would imply different population genetic forces, but since the differences in divergence times and distribution of taxa in the two subgroups / clusters being compared are not the same, the pi metric may simply reflect those historical contingencies. The metric pi is sensitive to the frequencies of alleles in a population, so more closely related mtDNAs will reduce pi, or more deeper-diverging mtDNAs will inflate it. So, some comment should be made to this effect to provide context for the larger pi values in the melanogaster subgroup. Notably, that group has three very diverged mtDNAs within the D. simulans clade that may inflate pi, which is not seen in the buzzatii sample. Clarify that the Dmel vs. Dbuzz difference is probably due to depth of divergence, not something distinct about rates of evolution. The use of a pairwise metric is fine, but its interpretation is complicated by the difference between a taxon sample vs. an allele sample from a breeding population.

Author’s response: we included a sentence explaining that the higher overall divergence observed for the melanogaster subgroup may be due to the deep divergence between some of the species pairs it includes, eg. D. melanogaster and D. yakuba, that diverged more than 10 MYA.

Lines 261-263: The higher dN/dS values for two genes in Dbuzz vs. Dmel data (Table2), may imply a shift in function. It could be noted that since the dS values (and pi as per points above) is higher in Dmel group, the true value of dS may be slightly underestimated due to multiple hits, leading to a slight overestimation of the dN/dS ratio in Dmel. This means that the elevated values of dN/dS for the ND and ATP genes in the Dbuzz group relative to the same Dmel genes are probably real.

Author’s response: We added a sentence (line 276) addressing the point raised by reviewer 2.

Submitted filename: Response_to_Reviewers.docx

Click here for additional data file.

2 Oct 2019

What does mitogenomics tell us about the evolutionary history of the Drosophila buzzatii cluster (repleta group)?

PONE-D-19-20340R1

Dear Dr. Moreyra,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

William J. Etges

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

28 Oct 2019

PONE-D-19-20340R1

What does mitogenomics tell us about the evolutionary history of the Drosophila buzzatii cluster (repleta group)?

Dear Dr. Moreyra:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. William J. Etges

Academic Editor

PLOS ONE