Bactrocera oleae ( Diptera : Tephritidae ) in Iran : An invasion from the Middle West

Despite an age-old tradition of olive growing and its geographical location, Iran was apparently free of the olive fly, Bactrocera oleae (Rossi, 1790) (Diptera: Tephritidae), the major worldwide olive tree pest, until the last decade. However, this situation has changed radically since B. oleae was first reported in 2004, and the olive fly is now Iran’s most important pest of olive trees, and a very serious threat to olive production in the Northern and more humid parts of the country. Genetic analyses of mitochondrial markers of B. oleae collected in the traditional olive growing area in Northwestern Iran were used to determine the origin of this pest. The results indicate it was most likely imported on olive trees from the Central Mediterranean area.


INTRODUCTION
The lands that presently integrate modern Iran have a millennia-long history of close ties with the Mediterranean.Hallmark Mediterranean agricultural crops like grapes and olives have a long history in the region, potentially as far back as the ancient civilizations of the Fertile Crescent (Djamali et al., 2009;Al-Ameri et al., 2011).Consistent with this view, a recent study found high levels of genetic diversity (and a clear divergence from other Mediterranean varieties) among Iranian olive trees (Hosseini-Mazinani et al., 2014), indicating a large population in the past and a distinct introduction into Iran.However, over the past two centuries, commercial olive tree cultivation has been largely restricted to the Northwest of the country, in the area of the Elburz Mountains (Hosseini-Mazinani et al., 2014).Because of this history and the country's geographic position it is surprising that the olive fly (Bactrocera oleae), the most important pest of olive trees worldwide, was apparently absent from Iran until very recently.Indeed, the native, Old World, range of B. oleae encompasses countries immediately to the East (Pakistan) and West (Turkey) of Iran.Unfortunately for Iranian growers, the situation has changed radically since B. oleae was first reported in 2004 (Jafari & Rezaei, 2004), with its status changing from a quarantine pest to the most important pest of olive trees (Noori, 2011).
There are several possibilities for the origins of the Iranian olive fly populations.First, given the presence of the pest in neighbouring countries, it could have diffused spontaneously from the East or the West.This could have been facilitated by the large increase in the area of olive orchards in Iran in recent times (Ghamary et al., 2010).On the other hand, this recent increase in the area covered with olive trees could have facilitated a dramatic increase in the size of residual autochthonous populations present in isolated patches of wild, feral, or abandoned olive trees.Finally, given that the enormous increase in the area planted with olive trees implied importing a large numbers of trees, an accidental human-mediated introduction was also a possibility.To distinguish among these hypotheses, we analyzed samples of olive flies collected in the most important olive growing region in Iran using mitochondrial DNA (mtDNA) markers capable of distinguishing between populations of different geographic origins (Nardi et al., 2005(Nardi et al., , 2010;;van Asch et al., 2012van Asch et al., , 2015;;Dogaç et al., 2013;Matallanas et al., 2013).

Olive fly samples
Olives were picked at eleven different locations (Table 1), seven in Iran, two in Italy, one in Tunisia and one in Palestine, and stored in plastic bags at room temperature.Emerging larvae, pupae and adults of B. oleae were collected periodically.Larvae were allowed to pupate in plastic containers, whereas adults and most pupae were stored at -20°C in 70% ethanol until DNA extraction.Some pupae were kept in plastic containers at room temperature until adults emerged.KT369102 + KT369105, respectively, while S1b3 haplotypes Iran01 and Iran03 can be recovered from the corresponding complete mtDNA sequences: KR677101 and KR677106.

Sequence and statistical analyses
Electropherograms were visualized using Chromas 1.45 (Technelysium Pty Ltd) and CLC Main Workbench 6 (CLC Bio).Sequences were manually corrected, aligned, cropped to the maximum overlap for each segment and assembled with CLC Main Workbench 6 (CLC Bio).File conversions were performed using FaBox (villesen, 2007) and Format Converter v2.2.5 available at http://www.hiv.lanl.gov/content/sequence/FORMAT_CONvER-SION/form.html.AMOvA, pairwise F ST values and haplotype and nucleotide diversities were calculated using Arlequin software v3.5.1.3 (Excoffier et al., 2005) and Tamura-Nei genetic distances (Tamura & Nei, 1993).For analyses of molecular variance, different populations, defined according to their geographic locations, were considered within each of the regional groups.Median-joining phylogenetic networks were constructed using the Network software v4.613 available at www.fluxus-engineering. com, using default parameters in all calculations (Bandelt et al., 1999).Models of nucleotide sequence evolution were selected with jModeltest (Darriba et al., 2012) using the Akaike Information Criteria with correction.Phylogenetic analyses were implemented in the BEAST1.8software package (Lemey et al. 2009;

DNA extraction, amplification and sequencing
Total DNA was extracted from individual flies using a standard SDS/Proteinase K method and used in separate amplification reactions.PCR amplification of the complete mtDNAs was achieved in fourteen segments (Table 2).For comparing haplotypes and phylogenetic analyses, amplifications and sequencing were restricted to "ND1p" -a previously analyzed 574bp fragment of the ND1 gene contained in segment S10 (Nardi et al., 2005;Dogaç et al., 2013), and "S1b3" -a combination of two regions (positions 522-1109 and 2513-3670 in the reference sequence -Genbank accession AY210702), contained in segments S01 and S03, totalling 1746bp (11% of B. oleae's complete mtDNA).Partial mtDNA sequences including S1b3 were used in previous studies (van Asch et al., 2012Asch et al., , 2015)), which facilitated comparisons with published results.Except for segment S11, PCR reactions were carried out in 25 μL reactions containing 10 ng of genomic DNA, 75 mM Tris-HCl (pH 8.8), 20 mM (NH 4 ) 2 SO 4 , 0.01% (v/v) Tween 20 (Fermentas, vilnius, Lithuania), 1.5 mM MgCl 2 (Fermentas), 0.25 mM of each deoxy-NTP (Fermentas), 25 pmol of each primer (Metabion International, Steinkirchen, Bavaria, Germany and Macrogen, Seoul, Republic of Korea) and 2.5 U of Taq DNA polymerase (Fermentas or Nzytech, Lisboa, Portugal).For S11, the concentrations of MgCl 2 and dNTP were doubled.The cycling protocol for segment 5 was: 95°C for 5 min; 3 cycles of 95°C for 30 s, 61°C for 1 min and 72°C for 1 min 45 s; 3 cycles of 95°C for 30 s, 58°C for 1 min and 72°C for 1 min 45 s; 3 cycles of 95°C for 30 s, 55°C for 1 min and 72°C for 1 min 45 s; 38 cycles of 95°C for 30 s, 58°C for 1 min and 72°C for 1 min 45 s; and 72°C for 5 min.For all the other segments, it was: 95°C for 5 min; 40 cycles of 95°C for 30 s, 58°C for 1 min and 72°C for 1 min 30 s; and 72°C for 5 min.Following treatment with ExonucleaseI (Fermentas) and Shrimp Alkaline Phosphatase (Fermentas), PCR products were Sanger-sequenced by Macrogen Europe (Amsterdam, The Netherlands).Primers used in amplification or sequencing are listed in Table 2. Complete mtDNA sequences are deposited in Genbank, with accession numbers KR677101 to KR677107.Of the newly identified ND1p haplotypes, H60 can be recovered from the corresponding complete mtDNA sequence (KR677101), whereas H61 to H69 are deposited in Genbank with accession numbers KR677092 to KR677100.Of the novel S1b3 haplotypes, Iran02, Iran04 and Iran05 can be reconstituted by joining sequences with accession numbers KT369100 + KT369103, KT369101 + KT369104 and   (Drummond et al., 2006).A coalescent tree prior with exponential population growth and priors on overall mutation rate, population size and population growth rate were chosen on the basis of previous results (Nardi et al., 2010).To verify convergence, two independent runs of 20 million generations (with a burn-in of 1 million and sampling every 1,000) were performed for each analysis.Trees were summarized and annotated using TreeAnnotator v1.8 and drawn using FigTree 1.4 (http://tree.bio.ed.ac.uk/ software/figtree/) and GIMP2.

RESULTS
Preliminary surveys conducted in the fall of 2013 revealed that although olive flies occurred as far east as Gorgan, in the Golestan province, high levels of infestation were only recorded in the Northwestern provinces of Qazvin, zanjan and Gilan.We therefore started by determining the complete mtDNA sequences of six individuals collected at different locations in these three provinces and one collected in Gorgan (Table 1 and Fig. 1), and compar- ing them to the 21 available sequences (Nardi et al., 2010).Bayesian phylogenetic analysis of these sequences confirmed the existence of three major clades corresponding to the three main areas of the original (known) Old World B. oleae geographic distribution: Pakistan, Sub-Saharan Africa and the Mediterranean basin (Nardi et al., 2010).Furthermore, the results obtained are consistent with the presence of three subclades in the Mediterranean basin, herein designated as P, M and O, as previously suggested (van Asch et al., 2015).Interestingly, the Iranian sequences clustered with the Mediterranean ones, indicating that neighbouring Pakistan could be excluded as a major source of current Iranian olive fly populations (Fig. 2).Furthermore, all seven Iranian sequences were found to belong to the same Mediterranean mtDNA clade (clade M).Since previous studies indicate that clade M predominates in the Central Mediterranean and is rare or absent in the Levantine area (Nardi et al., 2010;van Asch et al., 2012van Asch et al., , 2015)), this supports the possibility of a Central Mediterranean origin for Iranian olive fly populations.
To address this hypothesis, we focused on ND1p, a fragment of mtDNA encompassing 574 bp of the ND1 gene previously analyzed in samples from different areas of the Mediterranean basin (Table 1).Sequences from an additional 41 Iranian individuals, 15 individuals from the Central Mediterranean (Italy and Tunisia) and 10 from the Levant (Palestine) were therefore obtained and compared to a total of 296 previously reported sequences from the Levant, the Western (Aegean) provinces of Turkey and the Central Mediterranean (see map in Fig. 1) (Nardi et al., 2005(Nardi et al., , 2012;;Dogaç et al., 2013).Phylogenetic analysis of ND1p haplotypes revealed that Iran and the Central Mediterranean have similar haplotype distributions, which are markedly different from the haplotype distribution in the Levant, whereas Western Turkey displays an "intermediate" distribution (Fig. 3).This analysis also lead to the identification of ten new ND1p haplotypes (H60-H69).
The degree of differentiation among these four regional groups was further evaluated by determining pairwise F ST values and analysis of molecular variance (AMOvA).Pairwise F ST values strongly support the existence of genetic differentiation between all pairs of populations, with the exception of the Iran-Central Mediterranean pair (Table 3).Similarly, AMOvA analysis revealed that the difference between groups is not significant for the Iran-Central Mediterranean pair.In contrast, it is relevant (~13%) for the Iran-Western Turkey pair, whereas in the Iran-Levant pair the difference between groups accounts for more than half of the genetic variation (56.1%) (Table 4).
The haplotype network (Fig. 3) also revealed the presence in both the Central Mediterranean and Iran, of "Western European" haplotypes (Dogaç et al., 2013), which were not recorded in the 288 individuals from Western Turkey and the Levant.As four of the five sequences reported from Portugal, at the Western end of the Mediterranean basin and two of the five available French sequences (Nardi et al., 2005) belong to this "Western European" group, we further explored the hypothesis of a contribution from the Western Mediterranean to the invasion of Iran.The analy-ses of Iranian flies were therefore extended to S1b3 (Table 1), a combination of two segments previously used to study population structure of B. oleae in Western Europe (van Asch et al., 2012(van Asch et al., , 2015)).A comparison of the 50 new Iranian S1b3 sequences with those from Italy, France and Iberia (Nardi et al., 2010;van Asch et al., 2012van Asch et al., , 2015) ) led to the identification of five new S1b3 haplotypes and, together with the ND1p and complete mtDNA analyses, a total of 36 new mtDNA sequence variants (Table 5).More importantly, it showed that the three main Mediterranean mtDNA lineages (P, M and O) are present in Iran at frequencies similar to those found in Italy, but markedly different from those in France or Iberia (Fig. 4 inset).An analogous picture is revealed by pairwise F ST values (0.161 and 0.505, p < 0.001 for Iran-France and Iran-Iberia, respectively, vs 0.031, p = 0.49 for Iran-Italy) and an AM-OvA (Table 4).Indeed, whereas in the Iran-Italy pair the contribution of intergroup differences to genetic variation is not significant, it is important for comparisons involving areas located further to the West, approaching 14% and 50% for the Iran-France and Iran-Iberia pairs, respectively.These results confirm the lack of genetic differentiation between the Iranian and Italian populations and a high level of differentiation between the Iranian population and those from Iberia and France, apparently ruling out the latter two regions as the source of Iranian B. oleae.However, a more complex picture emerges when one considers P lineage substructure (Table 6).Indeed, clade P3 haplotypes, which are present in ~20% of Italian flies, were not recorded in Iranian samples, in which they are "replaced" by clade P1 haplotypes, typical of Iberia but with only a marginal presence in Italy.This discrepancy raises the possibility of an Iberian contribution to Iranian olive fly populations, implying the occurrence of more than one introduction event.However, an alternative explanation is that a population bottleneck at the time of the introduction led to alterations in haplotype lineage distributions.To evaluate this possibility, haplotype and nucleotide diversity levels in Iranian and Italian samples were analyzed.As shown in Table 7, both diversity levels were found to be much lower in Iran

DISCUSSION AND CONCLUSIONS
In the present study, we used molecular genetics to determine the source of the olive fly infestation that has struck Iran in the past decade.We analyzed mtDNA, because it has general features that make it an ideal choice for this type of work and, more importantly, it has been used to differentiate among B. oleae populations from different geographic locations, even within the Mediterranean basin (Nardi et al., 2005(Nardi et al., , 2010;;van Asch et al., 2012van Asch et al., , 2015;;Dogaç et al., 2013;Matallanas et al., 2013).In a pioneer study, Nardi et al. (2005) chose to analyze the ND1p 574bp fragment in a large number of individuals (n = 93) sampled from a variety of geographical locations worldwide.Their data seemed to indicate that this region of mtDNA could be used to distinguish populations from the three main areas of B. oleae's native range -Pakistan, Sub-Saharan Africa and the Mediterranean basin, but not from different parts of the Mediterranean basin.However, later work (Dogaç et al., 2013) has shown it can be used to differentiate between populations from the Western (Aegean) and Levantine provinces of Turkey.In a separate study, complete mtDNA sequences from a smaller set of individuals (n = 21) were used (Nardi et al., 2010), which, along with a Bayesian phylogenetic analysis, identified five well-supported mtD-NA clades, one formed by the Pakistani sequences, one by those from sub-Saharan Africa, and three by those from the Mediterranean area.Importantly, they were also able to distinguish between at least two groups of Mediterranean olive fly populations: Eastern (Levantine) and Central/ Western.In our previous studies (van Asch et al., 2012Asch et al., , 2015)), we opted for an "intermediate" strategy, based on using fractions of mtDNA ranging from 11.5% to 24.1%, but always including S1b3.This enabled us, not only to "recover" the five previously identified mtDNA sequence clades, but also to reveal the existence of well-supported subdivisions in the Mediterranean clades (van Asch et al., 2015).In addition, it showed that good discrimination among B. oleae populations from the Central and Western Mediterranean basin could be achieved, even when the sequence analysis was restricted to the S1b3 fraction of mtDNA (van Asch et al., 2012(van Asch et al., , 2015)).Interestingly, some differentiation between Iberian and Central Mediterranean populations was also detected using a 1151bp portion of the COI gene (Matallanas et al., 2013).
Because complete mtDNAs provide the best discriminating power, we decided to start this study by determining complete sequences of seven individuals collected at different locations in Iran.Though limited, this is the largest set of complete mtDNA sequences reported for a single country and provided important information regarding the origin of Iranian olive flies.Bayesian analysis strongly supported the notion that all seven sequences belong to the Mediterranean lineage M, essentially ruling out the possibilities that endogenous, Eastern (Pakistani) or Levantine populations were the main sources of Iranian flies.Indeed, the Pakistani sequences form a separate clade, estimated to have diverged from all other sequences 0.5 to 1.5 million years ago, which is consistent with the notion that Pakistani olive flies form a separate taxonomic entity (Silvestri, 1916;Nardi et al., 2010).On the other hand, recent genetic data indicate high levels of genetic diversity among Iranian olive trees and that the majority of the countrie's varieties cluster separately from Mediterranean cultivars (Hosseini-Mazinani et al., 2014).In light of these findings, one would expect that an endogenous population of olive flies would have a more diverse genetic structure than we found, with a clear separation, or at least the occurrence of mtDNAs occupying more basal positions on the phylogenetic tree, rather than the full integration into the Mediterranean clade observed.As for the Levant, considering its geographic position as the closest olive producing region West of Iran and the historical connections between the two areas, it was considered to be the main a priori candidate for the source of the olive fly in Iran.However, data from previous studies (Nardi et al., 2010;van Asch et al., 2015) indicate a strong predominance of the O lineage among Levantine populations, in contrast to the preponderance of M lineage in Iran.
The downside of opting for complete mtDNA sequences is that it usually (as was the case here) implies the use of a limited number of samples, therefore increasing the risk of sample bias (see below).Thus, more accurate genetic characterization of Iranian olive flies and comparisons with populations from other geographic areas were required, which involved carrying out mtDNA analyses of a larger number of individuals.Partial mtDNA sequences were therefore obtained from a total of 50 individuals from Iran, 15 from Italy and Tunisia and 10 from Palestine (to increase possible source population sample sizes from the Central Mediterranean and Southern Levant, respectively).Comparisons of populations from Iran, the Central Mediterranean, Western Turkey and the Levant using ND1p yielded not only 10 new haplotypes but, more importantly, provided evidence of genetic differentiation between all  et al., 2003).
In addition, haplotype network analysis highlighted the presence in Iran of sequences belonging to a group previously referred to as "Western European" (Dogaç et al., 2013), due to its prominent presence in both Portugal and France (Nardi et al., 2005).Since haplotypes from this group were not detected among a total of 163 sequences from Western Turkey or 125 sequences from the Levant, this result virtually excludes an Eastern Mediterranean source of Iranian olive flies.On the other hand, it raised the question: from "how far West" did the Iranian B. oleae come?To address this question, a different fraction of mtDNA, capable of distinguishing B. oleae populations from Western Europe (van Asch et al., 2012(van Asch et al., , 2015)), was analyzed.Both pairwise F ST and AMOvA-based comparisons between the Iranian population group on the one hand and those from Iberia, France and Italy on the other indicated that the Iranian group can be differentiated from those of Iberia and France, but not from the Italian group.This suggested that the Western Mediterranean region could be ruled out as a major source of Iranian flies, a hypothesis also supported by Bayesian analysis.Interestingly, Bayesian analysis of S1b3 sequences revealed that sequences belonging to clades P and O are also present in Iran, suggesting that complete mtDNA sequencing was able to detect only clade M due to an ascertainment bias.Nevertheless, the broader S1b3 analysis added further support for excluding Pakistan or the Levant as sources of Iranian olive flies.Indeed, Iranian sequences do not form separate clusters or tend to occupy basal positions in the phylogenetic tree (e.g.all of the Iranian M haplotypes are in the "core" Mn clade), none of them belong to the Pakistan clade and only one is in the Mediterranean clade O, which is dominant in the Levant (Nardi et al., 2010(Nardi et al., , 2015)).Instead, they are distributed among the three main Mediterranean clades, with frequencies similar to those found in Italy (see Fig. 3 inset), which is consistent with the hypothesis of a "Midwestern" origin for Iranian B. oleae.However, a discrepancy was observed between the distributions of sub lineages of Iranian and Italian clade P haplotypes, with all the Iranian lineages belonging to the P1 sub lineage, which is predominant in Iberia.This suggests that the possibility of more than one introduction, from different geographic origins, cannot be ruled out.However, this discrepancy could also be explained by contingency.Indeed, genetic diversity was much higher in Italy than in Iran, suggesting a population bottleneck at the time of the introduction that could have altered the distribution of P sublineages.
Regardless of whether one or multiple sources contributed to the olive fly populations currently infesting Iranian olive trees, the current work provides strong evidence that: (a) Central Mediterranean populations were the main source and, consequently, (b) that this was the result of human activity.This is the second case of human-mediated introduction of the olive fly in the last decade, underscor-ing the need for extremely tight control measures to prevent further accidental range expansions, particularly in the face of widespread international projects to expand olive culture, which will involve the transport of very large numbers of trees.In this regard, it is important to mention that Iranian farmers bitterly complained to us that they did not "have" olive flies until the new plantations were established -likely involving imported trees.This anecdotal evidence is supported by the data presented in this study, and the farmers' desperation is justified by the scale of the attacks we observed in loco.The recent introduction of the pest and the absence of natural enemies (similar to California) could have catastrophic consequences for traditional Iranian growers and even jeopardize the country's attempt to significantly increase olive oil production.A major effort, potentially involving the introduction of natural enemies, will be necessary to avoid such dire consequences and control the spread of olive fly in Iran.Some of the country's natural features, such as long periods of low rainfall and large areas of forest steppe and semi-desert might help to prevent B. oleae from spreading further, and even to eradicate it from some areas.However, given the milder climate in the traditional olive growing area of Iran, it seems likely that the olive fly will not be eliminated from there and that its control and management will have to become as much a part of olive tree agriculture as it is throughout the rest of the Mediterranean World.

Fig. 1 .
Fig. 1.Map showing the areas where samples of B. oleae were collected for determining the mtDNA sequences compared in this study.Regions that are genetically different are depicted in different shades of grey.The proposed introduction of Iranian populations from the Central Mediterranean is indicated by the arrow connecting the two regions, which are depicted in the same shade of grey.

Fig. 3 .Fig. 4 .
Fig. 3. Median-joining phylogenetic network comparing B. oleae ND1p haplotype distributions in Iran (green), the Levant (purple), Western Turkey (blue-grey) and the Central Mediterranean (dark blue).Circles represent haplotypes, black circles represent unobserved intermediate haplotypes, and the areas of the circles are proportional to the frequencies of the corresponding haplotypes.Haplotypes are numbered following Dogaç et al. (2013).The dashed oval indicates "Western European" haplotypes.

TaBle 1 .
Locations where Bactrocera oleae were sampled and the number of individuals used in each partial mtDNA analysis.

TaBle 4 .
Analyses of molecular variation (AMOVA) in the ND1p and S1b3 fractions of mtDNA used to compare Iranian olive fly populations and those belonging to the indicated groups.WTurkey -Western Turkey; CMed -Central Mediterranean.values represent the percentage of variance attributable to each source.ND1 S01b03 Iran vs Levant Iran vs WTurkey Iran vs CMed Iran vs Italy Iran vs France Iran vs Iberia

TaBle 7 .
(Rice, 2000;Riceiversity in the S1b3 fraction of B. oleae collected from Iran and six Italian locations.n-number of individuals, h -number of haplotypes, H -haplotype diversity, and π -nucleotide diversity.Central Mediterranean source of the Iranian olive flies, and a human-mediated invasion as in the case of North America(Rice, 2000;Rice