Introduction
Bluetongue (BT) and African horse sickness (AHS) are two non-contagious, arthropod-borne viral diseases affecting ruminants and horses, respectively. Bluetongue disease is caused by the bluetongue virus (BTV), the prototype of the Orbivirus genus within the Reoviridae family, which also includes AHSV (Holmes et al., Reference Holmes, Boccardo, Estes, Furuichi, Murphy, Fauquet, Bishop and Ghabrial1995). Both viruses are transmitted by species of biting midges belonging to the Culicoides genus (Diptera: Ceratopogonidae) (Du Toit, Reference Du Toit1944; Mellor & Pitzolis, Reference Mellor and Pitzolis1979; Mellor et al., Reference Mellor, Boned, Hamblin and Graham1990) and are maintained naturally through a series of alternative cycles of replication between Culicoides vectors and susceptible hosts (Takamatsu et al., Reference Takamatsu, Mellor, Mertens, Kirkham, Burroughs and Parkhouse2003).
In the Mediterranean basin and in southern Europe, the main vector for both viruses is Culicoides imicola (Mellor, Reference Mellor1996) although other Paleartic Culicoides species, including C. obsoletus, C. scoticus, C dewulfi and C. chiopterus, can also serve as main BTV vectors (De Liberato et al., Reference De Liberato, Scavia, Lorenzetti, Scaramozzino, Amaddeo, Cardeti, Scicluna, Ferrari and Autorino2005; Mehlhorn et al., Reference Mehlhorn, Walldorf, Klimpel, Jahn, Jaeger, Eschweiler, Hoffmann and Beer2007; Meiswinkel et al., Reference Meiswinkel, van Rijn, Leijs and Goffredo2007, Reference Meiswinkel, Baldet, de Deken, Takken, Delecolle and Mellor2008; Calvete et al., Reference Calvete, Calvo, Calavia, Miranda, Borrás, Estrada, Lucientes, Mañuz and Romero2008). In northern Europe, species of Culicoides obsoletus and scoticus, and the Avaritia subgenus (Culicoides dewulfi and chiopterus), are thought to be involved in bluetongue pothology.
In Europe, several BTV episodes during the past few decades have been reposted, with the first of these occurring between 1956 and 1960 in the southwestern region of the Iberian Peninsula (Manso-Ribeiro et al., Reference Manso-Ribeiro, Rosa-Azevedo, Noroña, Blanco-Forte, Grave-Pereira and Vasco-Fernandez1957). Since then, sporadic occurrences have been documented across Mediterranean Europe (Mellor & Pitzolis, Reference Mellor and Pitzolis1979; Jennings et al., Reference Jennings, Boorman and Ergün1983), with a large BT epizootic in 1998 expanding northward across Europe. During this epidemic, several BTV serotypes were postulated to affect many European countries, including several that had not previously been affected by the virus (Mellor & Wittmann, Reference Mellor and Wittmann2002; Mehlhorn et al., Reference Mehlhorn, Walldorf, Klimpel, Jahn, Jaeger, Eschweiler, Hoffmann and Beer2007).
Between 1956 and 1990, one BT and two AHS epizootics broke out in Spain, concentrating mainly in Andalusia (southern Spain) (Mellor & Pitzolis, Reference Mellor and Pitzolis1979; Jennings et al., Reference Jennings, Boorman and Ergün1983; Ortega et al., Reference Ortega, Holbrook and Lloyd1999); and, more recently, a BT epizootic in 2000 occurred in the Balearic Islands of Majorca and Menorca (Office International des Epizooties, 2000). Subsequent trapping revealed that C. imicola was the most likely vector involved, with Culicoides obsoletus and scoticus also abundantly present on these two islands (Miranda et al., Reference Miranda, Borras, Rincon and Alemany2003). C. imicola was first detected in 2000 on the Mediterranean Islands of Sardinia, Sicily and Corsica, located east of the Balearic Islands, as well as in the Italian mainland regions of Calabria and Tuscany (Wittmann et al., Reference Wittmann, Mellor and Baylis2001). BTV-4 outbreaks were reported in Menorca in 2003 and in Majorca by 2004, which then subsequently spread during the following year to other southern and mainland municipalities. In 2007, another BTV-1 outbreak emerged on the southern coast of the Iberian Peninsula, most likely transported by a viraemic animal or infected Culicoides specimen from northern Africa. The distribution of BTV-1 spread was similar to that of C. imicola but also spread to more northern areas where this species is absent (http://www.oie.int).
This northeastward expansion of the BT virus across Europe appears to be partially linked to the recent trend of C. imicola across Mediterranean Europe. Extensive surveys estimating current C. imicola distributions have been carried out in Portugal, Italy and Spain (Calistri et al., Reference Calistri, Goffredo, Caporale and Meiswinkel2003; Capela et al., Reference Capela, Purse, Pena, Wittmann, Margarita, Capela, Romao, Mellor and Baylis2003; Calvete et al., Reference Calvete, Miranda, Estrada, Borras, Sarto, Monteys, Collantes, Garcia-de-Francisco, Moreno and Lucientes2006). Sarto I Monteys et al. (Reference Sarto I Monteys, Ventura, Pages, Aranda and Escosa2005) have suggested that C. imicola is currently expanding into the Catalonian region (northeastern Spain). Furthermore, Calvete et al. (Reference Calvete, Miranda, Estrada, Borras, Sarto, Monteys, Collantes, Garcia-de-Francisco, Moreno and Lucientes2006) detected a C. imicola presence throughout Spain, but the considerable place-to-place and year-to-year variation of C. imicola populations (Rawlings et al., Reference Rawlings, Pro, Peña, Ortega and Capela1997), in addition to the lack of prior extensive surveys, does not allow for an accurate assessment of whether this species is indeed expanding into the remaining areas of Spain.
Understanding population genetic structures of insect vectors has important implications for the spatial scales over which vector-borne diseases spread (Tabachnick & Black, Reference Tabachnick and Black1995). Dipteran disease vectors differ enormously in the degree of matrilineal subdivision. C. imicola carried by midges from Portugal, Rhodes, Israel and South Africa belong to the same phylogenetic clade of mitochondrial cytochrome c oxidase subunit I (COI) and are phylogenetically distinct from the four other species of the C. imicola species-complex (Dallas et al., Reference Dallas, Cruickshank, Linton, Nolan, Patakakis, Braverman, Capela, Capela, Pena, Meiswinkel, Ortega, Baylis, Mellor and Mordue2003). This suggests that the South African C. imicola is phylogenetically distinct (Linton et al., Reference Linton, Mordue, Cruickshank, Meiswinkel, Mellor and Dallas2002) and that the level of the C. imicola matrilineal subdivision between the eastern and western ends of the Mediterranean basin is remarkably high, and is comparable to the geographical structure of BTV serotypes observed in outbreaks around the Mediterranean between 1998–2001 (Dallas et al., Reference Dallas, Cruickshank, Linton, Nolan, Patakakis, Braverman, Capela, Capela, Pena, Meiswinkel, Ortega, Baylis, Mellor and Mordue2003). Recent phylogenetic studies of the Culicoides species in France and Italy have been performed based on nuclear ITS1-rDNA and ITS2-rDNA sequences (Perrin et al., Reference Perrin, Cetre-Sossah, Mathieu, Baldet, Delecolle and Albina2006; Gomulski et al., Reference Gomulski, Meiswinkel, Delecolle, Goffredo and Gasperi2005, Reference Gomulski, Meiswinkel, Delecolle, Goffredo and Gasperi2006; Nolan et al., Reference Nolan, Carpenter, Barber, Mellor, Dallas, Mordue, Mordue and Piertney2007).
The aim of the present study was to characterize the C. imicola population in Spain and its relationship with other C. imicola populations in other areas. By comparing differences in mitochondrial cytochrome c oxidase subunit I (COI) gene, we suggest a hypothesis regarding current C. imicola expansion in Spain.
Materials and methods
Insect samples
Trapping was carried out on farms using a 4 W ultraviolet light trap fitted with a suction fan (Miniature Blacklight Model 1212, John Hock Company, Gainesville, FL, USA). Each trap was placed outside selected farms within 30 m of livestock and 1.7–2.5 m from ground level. All collected insects were transported to the laboratory and preserved in 70% ethanol. Upon examination, species of Culicoides were separated from other insects, and identification of C. imicola was achieved based on previously described wing patterns (Rawlings, Reference Rawlings1996). Between summer of 2005 and summer of 2006, 215 insects were randomly selected from 58 sites in Spain (fig. 1, table 1). Insects from regions of the most recent BT outbreaks and where C. imicola had been previously found were sampled along with samples from sites where C. imicola had been absent (sites 47–58). In places with four insects or less, this number was the maximum number of insects trapped. Sample sites in Spain were approximately 6–1100 km apart, with the pairwise distance between the centre of the sampling area approximately 290 km.
Fig. 1. Location of farms sampled under the Spanish National Surveillance Programme for the detection of Culicoides imicola in the present study. Dark grey: provinces in which C. imicola had been previously reported; light grey: provinces in which C. imicola was reported for the first time. Population abbreviations are as follows: Mijas (01); Tarifa (02); Navaluenga (03); Tudela de Duero (04); Fuente del Fresno (05); Ituero de Azaba (06); Vera (07); Moguer (08); Badajoz (09); Cedilla (10); Begijar (11); Valera de Abajo (12); Ceuta (13); Mallorca (14); Sant Lluis (15); Orihuela (16); Villar del Infantado (17); Herrera del Duque (18); Serrejón (19); Huerta de Valdecarabanos (20); Arévalo (21); Banobarez (22); Mula (23); Montemayor (24); Hellín (25); Cabeza de Buey (26); Fuente del Maestre (27); La Puebla de Alcocer (28); Navalmoral de la Mata (29); Arcos de la Frontera (30); Alcubillas (31); Salobrena (32); Lora del Río (33); Calzada de Oropesa (34); Arenas de San Pedro (35); Candeleda (36); Ibarzos (37); Guadalajara (38); Aldea del Fresno (39); Caldes de Malavella (40); Piera (41); Lorca (42); Molina de Aragón (43); Albalate de Cinca (44); Fraga (45); Hecho (46); Castejón (47); Gándara (48); Estivada (49); Fabara (50); Torremocha de Ayllón (51); Tejada (52); Asturias (53); Zaragoza (54); Montanana (55); San Blas (56); Villarroya (57); and Iglesia (58).
Table 1. Haplotype sequences and distribution of Spanish Culicoides imicola populations. IMICOI indicate the COI C. imicola haplotype. IMICOI 24-IMICOI 34 represent new haplotypes described in this study. Vertical numbers indicate the SNP position relative to the DQ868882 reference sequence. Population abbreviations are as follows: Mijas (01); Tarifa (02); Navaluenga (03); Tudela de Duero (04); Fuente del Fresno (05); Ituero de Azaba (06); Vera (07); Moguer (08); Badajoz (09); Cedilla (10); Begijar (11); Valera de Abajo (12); Ceuta (13); Mallorca (14); Sant Lluis (15); Orihuela (16); Villar del Infantado (17); Herrera del Duque (18); Serrejón (19); Huerta de Valdecarabanos (20); Arévalo (21); Banobarez (22); Mula (23); Montemayor (24); Hellín (25); Cabeza de Buey (26); Fuente del Maestre (27); La Puebla de Alcocer (28); Navalmoral de la Mata (29); Arcos de la Frontera (30); Alcubillas (31); Salobrena (32); Lora del Río (33); Calzada de Oropesa (34); Arenas de San Pedro (35); Candeleda (36); Ibarzos (37); Guadalajara (38); Aldea del Fresno (39); Caldes de Malavella (40); Piera (41); Lorca (42); Molina de Aragón (43); Albalate de Cinca (44); Fraga (45); Hecho (46); Castejón (47); Gándara (48); Estivada (49); Fabara (50); Torremocha de Ayllón (51); Tejada (52); Asturias (53); Zaragoza (54); Montanana (55); San Blas (56); Villarroya (57); and Iglesia (58).
DNA extraction, polymerase chain reaction (PCR) amplification and DNA sequencing
Total DNA was extracted from single midges using a commercial kit (Nucleo Spin Tissue, Mancherey-Nagel, Düren, Germany). The mitochondrial COI gene was amplified using the following primers: 5′-GGAGGATTTGGAAATTGATTAGT-3′ and 5′- CAGGTAAAATTAAAATATAAACTTCTGG-3′ (Dallas et al., Reference Dallas, Cruickshank, Linton, Nolan, Patakakis, Braverman, Capela, Capela, Pena, Meiswinkel, Ortega, Baylis, Mellor and Mordue2003). PCRs were carried out in a final volume of 25 μl containing 50 ng of total DNA, 5 pmol of each primer, 200 μm dNTPs, 2 mm MgCl2, 50 mm KCl, 10 mm Tris-HCl, 0.1% Triton X-100 and 0.5 U Taq polymerase (Taq polymerase, Biotools, Spain). Thirty cycles were performed with the following step-cycle program: strand denaturation at 94°C for 30 s, primer annealing at 54°C for 30 s, and primer extension at 72°C for 30 s. A 523 bp fragment of the COI gene was amplified from DNA extracts of individual midges (1–10 midges from each site), and PCR products were directly sequenced using the primers described above. All amplified products were sequenced from both strands.
Analytical methodologies
Sequence traces of each PCR fragment were edited in BioEdit v7.1.0 (Isis Pharmaceuticals, Inc, USA) to retain the maximum length sequence reads in both forward and reverse directions. Data were aligned without gaps and inspected using MEGA 2.0 (Kumar et al., Reference Kumar, Tamura, Jakobsen and Nei2001) to discriminate between singletons (polymorphisms that appeared in only one insect) and parsimony informative sites. COI reading frames were detected by alignment with other Culicoides COI sequences (table 2). Indices of sequence variation, and haplotype structure and neutrality tests were then calculated using DnaSP 4.0 (Rozas et al., Reference Rozas, Sanchez-DelBarrio, Messeguer and Rozas2003). Tajima's D (Tajima, Reference Tajima1989) tests the difference of θπ and θw, both diversity estimators of the neutral parameter 4 Neμ. θπ is the average heterozygosity per nucleotide site obtained from the mean number of pairwise differences among sample sequences (Nei & Li, Reference Nei and Li1979), and θw is estimated from the expected nucleotide diversity and number of segregating sites (Watterson, Reference Watterson1975). If the population sample fit the infinite sites model, θπ and θw would have the same expected value and D=0. Under non-neutral evolution, since θw is calculated only using the number of segregating sites, it is more affected than θπ by low-frequency variants calculated using the frequency of each variant. Therefore, an excess of singletons (as expected in cases of recent selective sweeps or population expansion) will produce negative Tajima's D values (Braverman et al., Reference Braverman, Hudson, Kaplan, Langley and Stephan1995; Aris-Brosou & Excoffier, Reference Aris-Brosou and Excoffier1996). Fu's Fs tests whether the excess of rare alleles is significant in relation to the expected value under a mutation-drift model. The D* statistic is based on the difference between the number of singletons and total number of segregating positions, while the F* statistic is based on the difference between the number of singletons and θπ. Both the D* and F* statistics test whether mutations are selectively neutral. The mismatch distribution was calculated using Arlequin (http://anthropologie.unige.ch/arlequin/) and represents the number of differences that exist between every pairwise combination of haplotypes plotted as a function of frequency. Mismatch distribution and the raggedness index of the observed distribution defined by Harpending (Reference Harpending1994) were also calculated to estimate recent demographic (Rogers & Harpending, Reference Rogers and Harpending1992; Schneider & Excoffier, Reference Schneider and Excoffier1999; Excoffier, Reference Excoffier2004) and spatial expansions (Ray et al., Reference Ray, Currat and Excoffier2003; Excoffier, Reference Excoffier2004). Significance levels of all test statistics were evaluated by randomization tests as described in Arlequin. Median networks using Network version 4.1.1.2 (http://www.fluxusengineering.com) were constructed to investigate the relationship between haplotypes. A minimum spanning network was built based on the absolute number of nucleotide differences. Data were analysed at two different hierarchical levels to test the relationship between Spanish C. imicola (entire dataset isolated in this study) and worldwide C. imicola haplotypes (entire dataset isolated in this study plus those deposited in GenBank from other countries and locations) (table 2). Non-Spanish samples were Genbank sequences obtained from Italy, Morocco, Portugal, Rhodes, Israel and South Africa.
Table 2. World Culicoides imicola haplotypes (IMICOI) and their distribution within the Mediterranean basin and in Europe. In this study, Spain 2006 refers to the Culicoides imicola survey. GenBank references are as follows: Spain 2006 haplotypes: DQ868882, DQ868883, DQ868886-DQ868890, DQ868892, DQ868894, DQ871030; SpainFootnote 1: AJ867227, AJ867228; ItalyFootnote 1: AJ877224, AJ877225, AJ867232-AJ867234; MoroccoFootnote 1: AJ867226; PortugalFootnote 2: AJ549415-AJ549426; PortugalFootnote 1: AJ867229-AJ867231; RhodesFootnote 2: AJ549388-AJ549392; IsraelFootnote 1: AF078097; IsraelFootnote 2: AJ549393-AJ549414, AF080531, AF080532, AF080534, AF080535, AF078098-AF078100, South AfricaFootnote 3: AF069231-AF069233, AF069249.
1 Haplotypes corresponding to GenBank sequences.
2 Haplotypes corresponding to sequences isolated in Dallas et al. (Reference Dallas, Cruickshank, Linton, Nolan, Patakakis, Braverman, Capela, Capela, Pena, Meiswinkel, Ortega, Baylis, Mellor and Mordue2003).
3 Haplotypes corresponding to sequences isolated in Linton et al. (Reference Linton, Mordue, Cruickshank, Meiswinkel, Mellor and Dallas2002).
Results
COI sequences of Culicoides imicola
A 455 bp mitochondrial sequence was collected from each of 215 C. imicola samples from 58 sites in Spain (fig. 1, table 1). The Spanish C. Imicola population contained 19 haplotypes defined by 15 polymorphic sites (four polymorphic sites were singletons), including two non-synonymous substitutions at positions 69 and 136 and three nucleotide variants in one polymorphic site at position 125 (SNP (Single Nucleotide polymorphim) position is relative to the DQ868882 reference sequence). Considering worldwide C. imicola populations, 34 haplotypes were defined by 27 mutations, including five non-synonymous substitutions at positions 46, 69, 115, 136 and 253. A total of 11 previously undescribed haplotypes were found among the Spanish C. imicola populations (GenBank Accession Nos. DQ868882, DQ868883, DQ868886-DQ868890, DQ868892, DQ868894 and DQ871030). The 11 phylogenetically informative sites were used to calculate nucleotide diversity (π=0.00132±0.17×10−3) and haplotype diversity (Hd=0.428±0.043) within the Spanish population. When separately considering each sampling location, the maximum nucleotide diversity and haplotype diversity were located in Ceuta (13) (π=0.00405, Hd=0.888), Tarifa (02) (π=0.00299, Hd=0.607) and Mijas (01) (π=0.00352, Hd=0.900) and the minimum only where the IMICO 1 haplotype was found (table 1). The 21 phylogenetically informative sites (six were singletons) in the worldwide C. imicola populations were used to calculate nucleotide (π=0.00344±0.14×10−3) and haplotype diversity (Hd=0.596±0.035).
Population characteristics of Culicoides imicola within the Mediterranean basin
To investigate haplotype structure, differences were calculated between each pairwise combination of the 19 Spanish and 34 worldwide haplotypes. The distribution of these differences indicated one peak in Spainish C. imicola haplotypes while two peaks were found worldwide, indicating the presence of divergent haplotype groups (fig. 2). This finding was supported by the construction of median networks (fig. 3), containing one median vector (mv1), which was missing one haplotype between western and eastern Mediterranean populations. Mv1 represents possible unsampled sequences or extinct ancestral sequences. It is clear when studying only the torso of the network, that haplotypes were also grouped in western Mediterranean (IMICOI 1, 14, 17, 24, 28, 31, 33) and eastern Mediterranean (IMICOI 2, 5, 6, 9, 13), showing in the latter two groups: Central Mediterranean, with samples of Rhodes, Sicily and one of Israel, and eastern Mediterranean. The last median vector (mv2) is a haplotype that connects this group to the second one, formed from Israel haplotypes. Fifteen haplotypes radiated from the western torso of the haplotype network, creating a clear star-like phylogeny, and samples from South Africa were situated outside of the torso of the network.
Fig. 2. Frequency distribution of the number of sequence mismatches between pairwise combinations of Culicoides imicola mtDNA haplotypes in Spain and worldwide. Distinct peaks indicate two haplotypes groups (–•–, Spain; –▪–, world).
Fig. 3. Median joining network diagrams show the approximate location of the two major geographic regions from which Culicoides imicola haplotypes were isolated. The diagram illustrates both the relationship between haplotypes and their observed frequencies. Haplotype number corresponds to IMICOI haplotype. The diameter of H_6 (IMICOI 6) represents one animal. Mutated nucleotides position (relative to the reference sequence) is indicated with numbers.
Neutrality tests and mismatch distributions
Neutrality tests were applied to determine whether neutrality could be discarded and to verify that the Spanish population was indeed expanding. Segregating sites displayed an irregular frequency distribution pattern, which represented an excess of variants that are both non- and highly polymorphic. Fu's Fs was −21.57196 (P<0.001), indicating a clear deviation from the expected allele frequency spectrum in the case of an excess of rare alleles. Nevertheless, the observed excess of rare alleles does not imply that COI variation is high. In fact, nucleotide variation was markedly reduced in our study. Despite the low number of segregating sites due to the excess of low frequency variants, the number of alleles was actually higher than expected. Tajima's D statistic was −1.99 (significant, P<0.05), while the D* and F* statistics were −1.353 (not significant, P>0.10) and −1.921 (not significant, 0.10>P>0.05), respectively. Negative values of neutrality estimators indicate an excess of low frequency variants. Fu's Fs was statistically significant, providing evidence for an excess of rare alleles produced by recent mutations. Tajima's D and Fu's Fs could indicate a recent population expansion, as there is an excess of low-frequency haplotypes with all new and rare mutations. Since Fu and Li's and Tajima's tests combine non-synonymous and synonymous mutations, it is impossible to separate effects of heterogeneous (purifying selection) and homogeneous (population expansion) processes. An excess of low-frequency haplotypes is expected under purifying selection, as new mutations with deleterious effects are removed (Hahn et al., Reference Hahn, Rausher and Cunningham2002). Pairwise mismatch distributions were calculated to distinguish whether Tajima's and Fu's Fs parameters could indicate a population expansion. Pairwise mismatch distributions were unimodal for Spanish COI sequences within Spanish C. imicola populations, as expected after a sudden increase in population size. The distribution of mtDNA did not differ significantly from those expected from a sudden expansion model (sum of squared deviations under the least-squares approach=0.00038, P=0.650; Harpending's raggedness index=0.12, P=0.800) (fig. 4). In the same way, the haplotype distribution did not differ significantly from those expected from a spatial expansion model (sum of squared deviations under the least-squares approach=0.00038, P=0.800; Harpending's raggedness index=0.12, P=0.700) (fig. 4). Harpending's raggedness index was very low, indicating that the observed distribution was unimodal. This index takes larger values for multimodal distributions than unimodal or smoother distributions commonly found in a stationary population. Thus, analysis of mismatch distributions was consistent with the hypothesis that the sampled Spanish C. imicola population had undergone a recent expansion.
Fig. 4. Frequency distribution of sequence mismatches between pairwise combinations of observed and estimated Spainish Culicoides imicola mtDNA haplotypes under the sudden and spatial population expansion models (–◆–, observed; –•–, sudden expansion model; 
, spatial expansion model).
Discussion
In the present work, we characterized a population of C. imicola in Spain to detect putative expansion across Spain by comparing differences in the mitochondrial cytochrome c oxidase subunit I (COI) gene. In the Spanish C. imicola population, a large number of rare haplotypes separated by single nucleotide differences were revealed. The distribution patterns of DNA variation, genotypes and haplotypes represent the first evidence of rapid C. imicola population expansion in recent times.
The Mijas (01), Tarifa (02) and Ceuta (13) sampling sites exhibit the maximum haplotype and nucleotide diversities. Ceuta is a north African Spanish town on the border of Morocco, while the other two are located along the southern coast of the Iberian Peninsula near Morocco. This could indicate that these places are a point of incursion of C. imicola in the Iberian Peninsula, as haplotypes described in Morocco were also found in these locations. This is in accord with the points of the BTV-4 outbreak in 2004 (Gomez-Tejedor, Reference Gomez-Tejedor2004) and the BTV-1 outbreak in 2007.
Mismatch distributions of various haplotypes were consistent with results from Dallas et al. (Reference Dallas, Cruickshank, Linton, Nolan, Patakakis, Braverman, Capela, Capela, Pena, Meiswinkel, Ortega, Baylis, Mellor and Mordue2003), who demonstrated that COI haplotypes in C. imicola were remarkably high in the eastern and western Mediterranean basin, concordant with the geographical structure of BTV serotypes in outbreaks reported around the area. Two peaks were identified, indicating the presence of divergent haplotypes groups (fig. 2). One peak contained no mismatches and the second had four. The first peak corresponded to haplotypes described in western Mediterranean while the second one corresponded to eastern Mediterranean haplotypes, indicating two genetically distinct sources better than genetic drift associated with dependence on fragmented area. This was confirmed by median joining networks that showed one median vector, ‘mv1’, (one missing haplotype) between western and eastern populations. The median network (fig. 3) also demonstrated that the haplotypes observed in the western Mediterranean region were closely related with one another, creating a clear star-like phylogeny centred around the most frequent haplotype. This type of population structure is best explained by rapid population expansion after a bottleneck, when molecular markers may show quite uniform structures over large areas. In such cases, the haplotype and nucleotide diversities can be low, even if the population de facto has some breeding structure (Walker et al., Reference Walker, Moler, Buhlmann and Avise1998). Another explanation for the observed star-like population structure is a selective sweep, where one genotype is favoured by natural selection and has recently replaced all others (Rich et al., Reference Rich, Licht, Hudson and Ayala1998). Median joining networks confirmed previous data obtained by Dallas et al. (Reference Dallas, Cruickshank, Linton, Nolan, Patakakis, Braverman, Capela, Capela, Pena, Meiswinkel, Ortega, Baylis, Mellor and Mordue2003), who suggested that gene flow between western and eastern areas in the Mediterranean area was restricted to the central Mediterranean region because insects isolated in Sicilia grouped with western and eastern specimens. The postulated northern expansion was further supported by the statistically significantly negative Tajima's D and Fu's Fs values, as well as expected mismatch distributions of expanding populations (fig. 4). Furthermore, D* and F* were negative, which again reflects population growth with no background selection. Samples drawn from populations at demographic equilibrium usually exhibit a multimodal distribution, as it reflects the highly stochastic shape of gene trees, but the distribution becomes unimodal in populations that have passed through a recent demographic expansion (Slatkin & Hudson, Reference Slatkin and Hudson1991; Rogers & Harpending, Reference Rogers and Harpending1992) or range expansion with high migration rates between neighbouring populations (Ray et al., Reference Ray, Currat and Excoffier2003; Excoffier, Reference Excoffier2004). A population generally undergoes spatial expansion if the population range is initially restricted to a very small area, then increases over time and space. Based on simulations from Ray et al. (Reference Ray, Currat and Excoffier2003), a large spatial expansion in a mismatch distribution often produces the same signal as for a purely demographic expansion in a panmictic population, but only if neighbouring subpopulations exchange many migrants (50 or more). These results, in conjunction with the presence of one common haplotype in almost all sample populations, support the recent spatial and demographic population expansion of C. imicola into northern regions (fig. 1). When western Mediterranean and worldwide imicola populations were combined, significant C. imicola expansion was similary found. However, no significance was seen when only eastern haplotypes were considered, possibly due to small sample sizes for sites in the eastern Mediterranean region. We would like to note, however, that results from worldwide C. imicola analysis should be interpreted with caution, as non-Spanish insects samples were collected in different years and, thus, could produce a bias in genetic analysis.
Previous data, as well as new data from this study, demonstrate that C. imicola is indeed expanding northwards. Expansion was limited to the southwest region of the Iberian Peninsula prior to 2000 and, then, was observed in more northern areas, such as Catalonia, the Balearic Islands and, recently, in the Spanish Basque country (Miranda et al., Reference Miranda, Borras, Rincon and Alemany2003; Sarto I Monteys et al., Reference Sarto I Monteys, Ventura, Pages, Aranda and Escosa2005; Calvete et al., Reference Calvete, Miranda, Estrada, Borras, Sarto, Monteys, Collantes, Garcia-de-Francisco, Moreno and Lucientes2006, Reference Calvete, Calvo, Calavia, Miranda, Borrás, Estrada, Lucientes, Mañuz and Romero2008). Our results imply that C. imicola expansion was a rapid and recent phenomenon. Ecosystem modifications depend on complex relationships established between the environment and the biology and physiology of the vector, virus and vertebrate hosts, providing a possible reason for the rapid expansion of C. imicola. Future research is necessary to assess the potential for other Culicoides species to act as vectors, such as the Culicoides pulicaris, and the putative competency of C. imicola, C. obsoletus group (C. obsoletus and scoticus), as well as species of the Avaritia subgenus (C. dewulfi and chiopterus).
Acknowledgements
The authors thank the veterinary services of regional governments for their efforts in carrying out the samplings. The national survey was supported by the Spanish Ministry of Agriculture, Fisheries and Food.
