Mitochondrial DNA variation of Drosophila obscura (Diptera: Drosophilidae) across Europe

Drosophila obscura is a common fruit fl y that inhabits the temperate forests of Europe. While it is abundant in the north compared to other Drosophila, its density decreases southwards, where it is gradually replaced by other Drosophila species. This study describes variation in the mitochondrial Cyt b gene of D. obscura from several European populations. We observed a large number of haplotypes, together with the structuring of genetic variation. Genetic variation is higher in the west where O1 and related divergent haplotypes dominate. In the east, the O2 haplotype is most frequent, together with haplotypes that recently arose from it. In the central part of the species range, both O1 and O2 are equally present, along with many others. These data reveal signs of population expansions that probably happened earlier in the west, and more recently in the east. Though our conclusions are based on only one genetic marker, limiting the power of the analysis, the results imply either postglacial expansion from two unique sources or, more likely, eastwards stepping-stone expansion. This study adds important information on genetic variation and phylogeography to the obscure biology of D. obscura, a species that has the potential to become an interesting model in evolutionary biology and conservation genetics. * Current address: College of Life and Environmental Sciences, University of Exeter, Penryn Campus, Cornwall TR10 9FE, UK; email: M.A.Wallace@exeter.ac.uk INTRODUCTION Drosophila obscura, the nominate species of the obscura group, is a Palearctic fruit fl y whose distribution extends from Southern Europe to the central Fennoscandia, and from Western Europe to central Asia (Lakovaara & Saura, 1971; Brehm & Krimbas, 1991). In Europe, its relative abundance decreases from north to south (Brehm & Krimbas, 1991). In southern Finland, it is the most abundant Drosophila species (Lakovaara & Saura, 1971). Along with D. subobscura, D. obscura is the one of the two most common Drosophila species in Great Britain, although D. subobscura outnumbers it during most months (Shorrocks, 1975; Begon, 1978). In the Central Balkans, among obscura group species, it is much less abundant than D. subobscura (Stanić et al., 2002; Pavković-Lučić et al., 2012), but is readily collected at higher altitudes (long-term observaEur. J. Entomol. 119: 99–110, 2022 doi: 10.14411/eje.2022.011

tion of the authors). Among obscura group species it forms a monophyletic triad with the two most related species: D. ambigua and D. trisits (Bachmann & Sperlich, 1993;Gao et al., 2007). Although these three species live in sympatry, D. obscura is more common (Stamenković-Radak et al., 2003) especially at higher altitudes compared to the other two, which are rarely collected there (long-term observation of the authors).
Species from the obscura group, especially American D. pseudoobscura and European D. subobscura, have long served as models in evolutionary biology (Anderson et al., 1975;Krimbas, 1993;Powell, 1997;Schaeffer et al., 2003;Balanyá et al., 2006;Savić Veselinović et al., 2019). The genetic variation of natural populations of D. subobscura has been well described using phenotypic, chromosomal, and molecular markers (Krimbas, 1993;Pascual Mitochondrial DNA variation of Drosophila obscura (Diptera: Drosophilidae) across Europe INTRODUCTION Drosophila obscura, the nominate species of the obscura group, is a Palearctic fruit fl y whose distribution extends from Southern Europe to the central Fennoscandia, and from Western Europe to central Asia (Lakovaara & Saura, 1971;Brehm & Krimbas, 1991). In Europe, its relative abundance decreases from north to south (Brehm & Krimbas, 1991). In southern Finland, it is the most abundant Drosophila species (Lakovaara & Saura, 1971). Along with D. subobscura, D. obscura is the one of the two most common Drosophila species in Great Britain, although D. subobscura outnumbers it during most months (Shorrocks, 1975;Begon, 1978). In the Central Balkans, among obscura group species, it is much less abundant than D. subobscura (Stanić et al., 2002;Pavković-Lučić et al., 2012), but is readily collected at higher altitudes (long-term observa-default parameters (Thompson et al., 1994). The fi nal analysis included 185 sequences.
A Median Joining network (Bandelt et al., 1999) was calculated and plotted in R v4.1.0 (R Development Core Team, 2018), using the pegas package (Paradis, 2010). Epsilon was set to 0, to put the fewest possible median haplotypes in the network. A Cyt b sequence of D. ambigua (IFS R42, collected at Mt. Rtanj in 2014) was used as an outgroup. The network was plotted to show only single alternative mutational steps between median haplotypes.
Bayesian inference of phylogeny implemented in Beast v2.6.2 (Bouckaert et al., 2019) was used to infer tree topology and then to construct a phylogenetic tree. Sequences of closely related species D. ambigua and D. tristis (NCBI accession number: EF216284.1) were used as outgroups to root the tree. We selected the best-fi t substitution model for these data using likelihood ratio tests and Akaike information criterion (Akaike, 1973), implemented in jModelTest v2.1. 10 (Posada, 2008). The best-fi t model was a Tamura-Nei model of nucleotide substitution with a significant proportion of invariable sites (I) and gamma (G) distributed among-site rate heterogeneity (TnR + I + G). The Markov chain et al., 2001;Fragata et al., 2010;Savić Veselinović et al., 2019). Particularly interesting is the pattern of mitochondrial (mtDNA) variation found in natural populations of D. subobscura, which has proven to be an excellent model for studying the selective forces that maintain sympatric mtDNA variation (Jelić et al., 2015;Savić Veselinović et al., 2019;Kurbalija Novičić et al., 2020).
In contrast to D. subobscura, limited data is available on genetic variation in natural populations of its sympatric counterpart D. obscura. Chromosomal inversion polymorphism has been studied in 13 isofemale strains (IFSs) collected across the European continent (Brehm & Krimbas, 1991). The sample size was limited, but the number of detected inversions would suggest that chromosomal variation in D. obscura is similar to that in D. subobscura. Genetic variation has also been assessed by enzyme loci in multiple populations from Fennoscandia (Lakovaara & Saura, 1971) where a large numbers of polymorphic loci were observed. However, so far, mtDNA variation has not been assessed in this species. There is also a lack of D. obscura population studies that cover wide geographic areas. This species is easily collected in the wild, easily bred in the laboratory, and is a promising model for studying historical and adaptive processes that have shaped the genetic variation of natural populations of the Palearctic.
In this paper, we describe variation in an 893 bp long sequence of the mitochondrial Cytochrome b (Cyt b) gene from several European populations of D. obscura. We record signifi cant genetic differentiation among different regions of Europe, as well as different levels of withinpopulation variation. We discuss our fi ndings in light of colonisation from glacial refugia, admixture, and more recent processes that could have infl uenced the observed pattern of variation.

MATERIALS AND METHODS
Samples of D. obscura were collected from four European countries: Serbia, Finland, Germany, and Scotland, UK (referred to as four populations in the further text), covering a wide geographic range (Fig. 1). Samples from Serbia cover several distinct localities. Table 1 contains information on the collection sites, date of collection, and the number of specimens.
DNA was extracted using a method that enriches genomic DNA with mtDNA (Martinez et al., 1992). For Serbian samples, extraction was done from the F 1 progeny of the females collected in the wild. Extraction was conducted separately for the progeny of different females. For the rest of the samples, DNA was extracted from wild-caught individuals that had been kept in ethanol, and the fi nal step of alkaline lysis was excluded to obtain enough DNA. An 893 bp fragment corresponding to the mitochondrial Cyt b gene was PCR amplifi ed and sequenced with the following primers: Cyt b-F 5'-TTAT GGTT GATT ATTA CGAA-3' and Cyt b-R 5'-CAAA ACAT ATGC TTAT TCAA-3' (Gao et al., 2007). The PCR cycling conditions consisted of an initial denaturation step at 94°C for 3 min, 35 cycles: at 94°C for 50 s, 51.5°C for 1 min, and 72°C for 1 min; with a fi nal extension at 72°C for 3 min. Amplifi ed products were purifi ed using the QIAquick PCR Purifi cation kit (QIAGEN, Hilden, Germany), and sequencing reactions were performed using both primers (Macrogen inc. Amsterdam, The Netherlands). The obtained sequences were aligned in BioEdit v.7.2.5 (Hall, 2011), using the ClustalW algorithm with  We then calculated nucleotide and haplotype diversity. Tajima's D (Tajima, 1989), Fu's Fs (Fu, 1997, Fu andLi's D-F (Fu &Li, 1993), andFay andWu's H (Fay &Wu, 2000) tests were used to test departure from mutation-drift equilibrium. We used several tests due to their differing statistical approaches, and power to infer demographic history (Ramos-Onsins & Rozas, 2002). This was important since populations differed in the number of analysed sequences. For Fu and Li's and Fay and Wu's tests, a Cyt b sequence from the most closely related species D. ambigua (IFS R42, collected at Mt. Rtanj in 2014) served as an outgroup. A McDonald-Kreitman test (McDonald & Kreitman, 1991) was implemented to compare the ratio of non-synonymous to synonymous change within and between species for D. obscura and D. ambigua. Changes in population size were examined by calculating the observed and expected pairwise differences (mismatch distribution) (Rogers & Harpending, 1992). The above parameters and tests were conducted using DNASP v.6.0 (Librado & Rozas, 2009). An extended Bayesian skyline plot (EBSP) (Heled & Drummond, 2008) was applied to additionally infer demographic history using BEAST2 (Bouckaert et al., 2019), with an appropriate substitution rate for D. melanogaster mtDNA (Haag-Liautard et al., 2008), and an assumption of 4 generations per year (Begon, 1978) which corresponds to the clock rate of 0.248. The length of the Markov chains was set to 10 9 for the EBSP, logging the parameters every 3000 iterations. Burn-in was set to discard 25% of the samples. Tracer (Rambaut et al., 2018) was used to assess the convergence of the chains.
An analysis of molecular variance (AMOVA) was implemented in order to partition variation between and within populations that represent the four geographical regions. An AMOVA was also conducted separately for the Serbian population that is comprised of four localities (with the exclusion of MS locality where only two individuals were available). Pairwise F ST indices were calculated using Arlequin v.3.5.1.2 (Excoffi er & Lischer, 2010). The Mantel test (Mantel, 1967) in Arlequin v.3.5.1.2 was used to test for a signifi cant isolation-by-distance correlation. The Mantel test compares F ST genetic distances with log e-transformed Euclidean spatial distances in kilometers. Another form of mantel test was employed, using two matrices of log e coordinates, one for longitude and one for latitude, to be compared to the F ST matrix.
All samples were tested for the presence of Wolbachia. A PCR assay using 16S rDNA Wolbachia-specifi c primers (O'Neill et al., 1992) was used with a slight modifi cation to the PCR conditions (García-Martínez et al., 1998). Drosophila tristis with known infection status served as a positive control (Erić et al., 2019). To exclude the possible presence of other maternally-transmitted bacteria (Hurst & Jiggins, 2005) we conducted microbiome sequencing. Two samples were made, each using 10 females originating from 10 different randomly chosen IFSs from Serbia. One sample included the O1 haplotype and 9 others that share its specifi c substitution on position 828. The other included the O2 hap-lotype and 9 others with its specifi c variant on position 828. At the time of the analysis, these lines had been kept in the laboratory for 5 generations. DNA was extracted from pooled individuals using a previously published protocol for microbiome sequencing in Drosophila (Kapun et al., 2020). The microbiome sequencing was performed by Fisabio (Valencia, Spain) and included Illumina 16S V3-V4 amplicon library preparation and MiSeq 300 bp paired-end sequencing. Primer sequences [forward: CCT ACG GGN GGC WGC AG, reverse: GAC TAC HVG GGT ATC TAA TCC (Klindworth et al., 2013)] were removed with BBDuk (https://jgi.doe.gov/data-and-tools/bbtools/) using kmer length 15 and allowing for 2 mismatches. The generated data was analysed in R v4.1.0 (R Development Core Team, 2018) using dada2 for error estimation, sequence denoising, merging, and chimera removal (Callahan et al., 2016) as described previously (Beribaka et al., 2021). In short, sequences were fi rst trimmed to 240 bp for forward and 210 bp for reverse reads, and all sequences containing more than 2 or 4 expected errors (for forward and reverse reads respectively) were discarded. After denoising, sequence pairs were merged using a minimum overlap of 20 bases without mismatches, and all sequences shorter than 400 or longer than 428 were discarded. Taxonomy assignment up to genus level was performed using the Silva v132 (https://www.arb-silva.de/documentation/release-132/) database with the RDP Naive Bayesian Classifi er algorithm (Wang et al., 2007) as implemented in dada2 with default parameters.

RESULTS
Among the 185 obtained sequences, we detected 72 different haplotypes for the Cyt b gene (NCBI accession numbers: MZ337620.1-MZ337804.1). There were 73 segregating sites, among which two had three alternating nucleotides. Among all segregating sites, 59 were synonymous, and 16 non-synonymous substitutions (Table S1). Moreover, 43 changes were observed only once, while 32 were observed in multiple individuals. One specifi c substitution (G to A on position 828 which changes valine to isoleucine) is of particular interest due to its presence in 73 individuals. The variant with valine is particularly frequent in Serbia (0.875), followed by Finland (0.750), and Germany (0.524), whilst it is rare in Scotland, UK (0.083). The lowest diversity was recorded in Serbia (π = 0.0146, hd = 0.546), while the highest was in Scotland (π = 0.0387, hd = 0.964, Table 1). The majority of haplotypes were singletons.
A Median Joining network of D. obscura Cyt b haplotypes is presented in Fig. 2. The O1 haplotype is the most common in Scotland, together with a large number of very divergent haplotypes that share its specifi c substitution on position 828. Also, Scottish haplotypes are more closely related to the outgroup species D. ambigua sequence (Fig.  2). The O2 haplotype is particularly frequent in Serbia, where most other haplotypes share its specifi c substitution on position 828 and differ in only one mutational step from O2. In addition to the O1 and O2 haplotypes, O3 is also observed in more than two sampled individuals. O3 shares a specifi c substitution at position 828 with the O1 haplotype. Except for O1, which was not observed in Finland, these three most frequent haplotypes are present in all analysed regions. Haplotype O2 is the most frequent overall, maybe due to the overrepresentation of samples from Serbia.
Apart from splitting outgroup species as different clades, the Bayesian tree ( Fig. S1) shows one sequence from Scotland (O65) as a monophyletic clade (PP = 1), and a few shallow nodes splitting sister haplotypes (PP > 0.7). After that, the phylogenetic relationships could not be resolved confi dently, as the posterior probability is very low for the rest of the nodes (PP < 50%).
The results of Tajima's D, Fu's Fs, Fu and Li's D-F, and Fay and Wu's H are presented in Table 2. Negative values were observed in all cases, the majority of which were statistically signifi cant. The McDonald-Kreitman test (Table  3) showed that the ratio of non-silent to silent variation was greater within species than between species. A statistically signifi cant departure from neutrality was observed for the Serbian population (SR) and the German (GF) population.
Graphs of mismatch distribution are presented in Fig. 3. All graphs show some raggedness for the observed function. A unimodal distribution with a peak at 0 observed pairwise differences is found for sequences from Serbia. A unimodal distribution with a peak at 2-3 observed pairwise nucleotide differences is found in Scotland. A bimodal distribution was observed for sequences from Finland, Germany, and total Europe. Visual inspection shows that only samples from Scotland show a good fi t to a population expansion scenario in contrast to constant population size.  EBSP estimates of the time of demographic expansion are given in Fig. 4. The analysis with all sequences included shows that population expansion started around 3000 years ago (Fig. 4a). When analysing the Serbian population the graph gives signs of recent expansion that started around 1200 years ago but with notable uncertainty (Fig.  4b). The EBSP plot for samples from Scotland estimates population expansion to have started at around 8000 years ago, but the confi dence interval gets wider as we go back in time (Fig. 4c). The German sample was very small and the EBSP function did not detect past population expan-sion, as the Markov chain did not converge even when we increased the number of iterations (Fig. 4d). The population from Finland, despite the small sample size, showed signs of a mild population expansion ranging from 3000 to 2000 years in the past, though with a wide confi dence interval (Fig. 4e).
The results of the AMOVA show that 13.24% of the total variation is present among populations (F ST = 0.1324, P < 0.001). Pairwise F ST values and the signifi cance of their difference from zero are presented in Table 4. The most differentiated populations are Serbia and Scotland, while  Signifi cant isolation by distance is not observed in our data set. When geographical distances are used in the Mantel test, a positive correlation was observed (r = 0.304757; P = 0.330740), which was not statistically signifi cant. If latitude is considered, the correlation is negative, and also non-signifi cant (r = -0.234862; P = 0.667790). The highest correlation was observed for longitude (r = 0.735025; P = 0.082810), which, though not statistically signifi cant, shows signs of possible population differentiation along the East-West cline.
All samples were negative for the presence of Wolbachia. Microbiome sequencing of available IFSs shows a complete absence of maternally transmitted bacteria (not a single read was recorded for Wolbachia, Spiroplasma, Microsporidia, or Rickettsia) that could infl uence mtDNA variation (Fig. 5). Lactobacillus together with Acetobacter comprise more than 99% of the microbiota in both samples while Ralstonia comprises less than 0.5% and all the other genera represent less than 0.1%.

DISCUSSION
In this study, we analyse nucleotide variation of the mitochondrial Cyt b gene in a widespread European Drosophila species, throughout its range. Our results show that D. obscura possesses a high level of mtDNA variation both within and between populations. The nucleotide diversity within populations observed for the D. . Sequences form a complex haplotype network with several star-shaped subnetworks, with a geographical structuring of genetic variation across Europe, mostly across the East-West axis. Interestingly, as O1 and O2 frequencies differ in the East and West, so do the frequencies of the other haplotypes that share their specifi c variants on the 828th nucleotide of the Cyt b gene. Variation is especially high in the western part of this species range. The haplotypes sharing the O1 variant on position 828 are more divergent from one another than the haplotypes sharing the O2 variant. Haplotypes sharing the O2 828 variant, in all but one case, differ from O2 by only one mutational step. Haplotypes O1 and O2 show similar frequencies in Central Europe. Although the number of individuals sampled from this area was restricted, and only one site from Germany was sampled, it is intuitive to conclude that the haplotype distribution in Central Europe is somewhere between eastern and western populations. The observed geographical structuring of mtDNA variation in D. obscura is especially interesting when compared to the sympatric D. subobscura, which shows geographic homogeneity in mtDNA variation across its entire range. During the past 40 years, studies of more than 30 populations of D. subobscura have shown that two main haplotypes are almost equally present in all populations and that there are less frequent population-specifi c haplotypes (Kurbalija Novičić et al., 2020). Cytochrome b is highly conserved and the sequences are very similar, so posterior probabilities were low for all but a few nodes. Bayesian inference did not reconstruct the Cyt b phylogeny with high confi dence, but we did draw a few conclusions from the analysis. The Bayesian tree shows that haplotypes that share a polymorphism at position 828 with the O1 haplotype, and that are frequent in the west, are indeed older and closer to Drosophila ambigua. The tree also shows that sequences that share the 828 O2 polymorphism are grouped together, but with a low posterior probability.
The observed excess of singletons and a large number of segregating sites are responsible for the negative values for the array of parameters that measure the departure of haplotype distribution from mutation-drift equilibrium. This result implies either population expansion or purifying selection. Positive values of neutrality indices from the McDonald-Kreitman test indicate an excess of non-silent polymorphism compared to divergence. This excess of amino acid polymorphism in the mtDNA, relative to divergence, is generally present in mice and humans (Nachman et al., 1994(Nachman et al., , 1996Templeton, 1996), and also in Drosophila (Kaneko et al., 1993;Ballard & Kreitman, 1994;Rand et al., 1994), particularly for the Cyt b gene (Ballard & Kreitman, 1994;Erić et al., 2019). The observed pattern of non-silent polymorphism can be interpreted (Rand & Kann, 1996, 1998 in the light of the nearly neutral model (Ohta, 1992a, b) that predicts the accumulation of mildly deleterious alleles that persist for short time within a population and do not contribute to divergence (Nachman, 1998;Weinreich & Rand, 2000;Meiklejohn et al., 2007). Older, slightly deleterious haplotypes cannot become progenitors of new lineages due to natural selection (Grant, 2015). This could be the case for the haplotype pattern observed in Scotland. Recent increases in effective population size can also generate artefactual evidence of positive selection if substitutions are slightly deleterious and if there is no selection upon synonymous codon use (Eyre-Walker, 2002). This may be the case for the O2 haplotype, and other haplotypes characteristic for Eastern Europe, whose network refl ects only recent expansion.
Mismatch distribution analysis shows unimodal distributions for the Serbian population, with low mismatch values. Unimodal distributions with high mismatch values are detected for Scotland. A unimodal distribution of mismatches, which refl ects demographic expansion, moves to higher values as mutations accumulate over time in a population (Grant, 2015), which implies earlier population expansion in the western range of the species.
The analysis of infection status with maternally transmitted microorganisms is important since they share the same mode of inheritance with mtDNA. They can often confound the inference of evolutionary history obtained by mtDNA markers, as variation is altered by selection acting on these microorganisms (Hurst & Jiggins, 2005). Importantly, all our samples were negative for Wolbachia. Additionally, we conducted microbiome sequencing on the Serbian samples, and maternally transmitted bacteria were also excluded (Spiroplasma, Microsporidia, and Rickettsia). Although we cannot exclude their presence in non-tested samples, the main haplotypes that are present Europe-wide were covered within our sample. While Lactobacillus and Acetobacter, which were detected in our samples, are common bacteria also found in IFSs of D. subobscura and D. melanogaster kept in the laboratory, D. obscura shows less diversity of bacterial genera compared to these two species (Beribaka et al., 2021).
What historical processes could have shaped the observed pattern of mtDNA variation in D. obscura? The different compositions of haplotypes in the east and the west could imply postglacial colonisation from at least two different sources with an admixture in central Europe (Fig. 6a). The Balkan Peninsula would be the hypothetical source of the O2 haplotype, while the source of the O1 haplotype could be the Iberian or Apennine peninsula. The Central European population perfectly matches this scenario, as it possesses eastern and western lineages in almost equal frequency. In addition, graphs of mismatch distribution show two peaks that correspond to the peaks of two distinct lineages. In that respect, postglacial colonisation of this species would follow the hedgehog's scenario (Hewitt, 2004). Although our EBSP intervals for population expansion are very wide, the results imply that the expansion probably happened earlier in the West than the East. Then the question arises as to why expansion after the last glacial maximum happened at different times from different glacial refugia, and why there is such a difference in the level of variation between different refugial sources?
One might imagine another scenario, of postglacial expansion from western peninsulas, or western cryptic refugia, and then gradual colonisation of the European continent by eastwards stepping-stone range expansion (Fig.  6b). In the East, the O2 haplotype might have increased its frequency to the detriment of other haplotypes due to the bottleneck effect, and subsequently generated an array of young singleton haplotypes in its recent expansion. In this scenario, postglacial colonisation of this species would be more similar to the brown bear's scenario (Hewitt, 2004). This hypothesis is supported by the earlier expansion of D. obscura in the west, approximately 8000 years ago, while the O2 haplotype network might refl ect only recent expansion after its colonisation of the highlands of the Balkan Peninsula. These expansions may also include those that generally happen yearly after winter or summer contractions, which are also observed in D. subobscura (Castro et al., 2010;Christie et al., 2010;Erić et al., 2019). Signs of ancient historical expansions may also have been masked over by annual contractions and expansions. It is also important to understand that the time of the expansion generated by EBSP used in this study is just a rough estimation based on a constant number of generations per year (Begon, 1978) which could have varied through time and across different regions.
We should also stress a specifi c aspect of this species' biology, which could be responsible for the observed pattern of structure, and level, of genetic variation. In the Balkan Peninsula, D. obscura is found in great numbers only in higher altitudes, while moving to the north it is frequently found in the lowlands too. The fact that it is the most common Drosophila species in southern Fennoscandia (Lakovaara & Saura, 1971) implies that it is adapted to colder climates. So, it is expected for southern D. obscura to show different population size dynamics compared to northern Europe. On the other hand, ecological conditions could be similar since vegetation and climate in higher altitudes in the south are similar to those from northern Europe. Populations from the south are more isolated and prone to genetic drift. Although populations can reach high numbers of individuals in the summer, they are restricted to small geographical areas and less prone to gene fl ow in the optimal part of the season that could restore variation.
To more accurately decipher the population history that has caused the observed pattern of variation in D. obscura mtDNA more sampling is needed, especially from the Apennine and the Iberian Peninsula, where D. obscura is more likely to be found in higher altitudes. Further analysis should also be conducted with additional genetic markers, both mitochondrial and nuclear (Brito & Edwards, 2009). This species' distribution in higher altitudes and latitudes also makes it a potential model for studying evolutionary change due to global warming. Although Drosophila fl ies can migrate easily, woodland habitats suitable for D. obscura cannot be formed at the pace dictated by global warming. Additionally, in some regions, there is not much space left for the fl ies to move, either to higher altitudes or northwards. Being very abundant, and easy to collect (particularly in northern Europe), designate and maintain in the laboratory, together with its rich mtDNA polymorphism, makes this species a promising model for these sorts of evolutionary studies.

REFERENCES
Haplotype