Medical and Veterinary Entomology (2014), doi: 10.1111/mve.12042

Phenotypic differentiation and phylogenetic signal of wing shape in western European biting midges, Culicoides spp., of the subgenus Avaritia ˜ OZ-MUN ˜ O Z 1 , S. T A L A V E R A 2 , S. C A R P E N T E R 3 , F. M U N 4 S. A. N I E L S E N , D. W E R N E R 5 and N. P A G E` S 2 1

Departament de Biologia Animal, de Biologia Vegetal i d’Ecologia, Facultat de Bioci`encies, Universitat Aut`onoma de Barcelona, Bellaterra (Cerdanyola del Vall`es), Barcelona, Spain, 2 Centre de Recerca en Sanitat Animal (CReSA), UAB-IRTA, Campus de la Universitat Aut`onoma de Barcelona, Bellaterra (Cerdanyola del Vall`es), Barcelona, Spain, 3 Vector-borne Viral Diseases Programme, The Pirbright Institute, Pirbright, Woking, Surrey, U.K., 4 Department of Environmental, Social and Spatial Change, Roskilde University, Roskilde, Denmark and 5 Leibniz-Centre for Agricultural Landscape Research, M¨uncheberg, Germany

Abstract. In the past decade biting midges of the subgenus Avaritia (Diptera: Ceratopogonidae) have been popular subjects of applied entomological studies in Europe owing to their implication as biological vectors in outbreaks of bluetongue and Schmallenberg viruses. This study uses a combination of cytochrome oxidase subunit I barcode sequencing and geometric morphometric analyses to investigate wing shape as a means to infer species identification within this subgenus. In addition the congruence of morphological data with different phylogenetic hypotheses is tested. Five different species of the subgenus Avaritia were considered in the study (C. obsoletus (Meigen); C. scoticus Kettle and Lawson; C. chiopterus (Meigen); C. dewulfi Goetghebuer and C. imicola (Kieffer)). The study demonstrated that over 90% of individuals could be separated correctly into species by their wing shape and that patterns of morphological differentiation derived from the geometric morphometric analyses were congruent with phylogenies generated from sequencing data. Morphological data produced are congruent with monophyly of the subgenus Avaritia and the exclusion of C. dewulfi from the group containing C. obsoletus, C. scoticus and C. chiopterus. The implications of these results and their importance in a wider context of integrating multiple data types to interpret both phylogeny and species characterization is discussed. Key words. Avaritia, COI, geometric morphometrics, phylogenetic signal, wing

morphology

Introduction The genus Culicoides Latreille (Diptera: Ceratopogonidae) currently contains more than 1300 extant species around the world (Borkent & Wirth, 1997; updated 2011). This genus has been extensively studied because it includes biological vectors of internationally important arboviruses affecting domestic livestock (Mellor et al., 2000). During the last decade, major

outbreaks of bluetongue and Schmallenberg viruses have had a significant impact on European livestock trade, both through direct disease losses and through movement restrictions to reduce spread (Saegerman et al., 2008; Hoffmann et al., 2012). In Europe, a large proportion of confirmed and implicated vectors of Culicoides-borne livestock arboviruses are placed within the subgenus Avaritia. These include the major afrotropical vector C. imicola and four different species commonly

Correspondence: Francesc Mu˜noz Mu˜noz, Departament de Biologia Animal, de Biologia Vegetal i d’Ecologia, Facultat de Bioci`encies, Universitat Aut`onoma de Barcelona, E-08193 Bellaterra (Cerdanyola del Vall`es), Barcelona, Spain. Tel.: +34 93 581 37 45; Fax: +34 93 581 13 21; E-mail: [email protected] © 2014 The Royal Entomological Society

1

2 F. Mu˜noz-Mu˜noz et al. referred to as the Obsoletus group or complex (C. obsoletus, C. scoticus, C. dewulfi and C. chiopterus; Saegerman et al., 2008). Some of these vector species are challenging to separate using classic morphological methods (Meiswinkel et al., 2004). As correct identification of disease vectors is essential in transmission risk assessment, recent efforts have centred upon developing molecular and morphological tools that allow a reliable discrimination between similar species of the Obsoletus group (Cˆetre-Sossah et al., 2004; Pag`es & Sarto i Monteys, 2005; Mathieu et al., 2007; Nolan et al., 2007; Augot et al., 2010; Nielsen & Kristensen, 2011). In spite of the advantages of employing morphometrics in cryptic species identification and the widespread use of this methodological approach in other biological disciplines, few studies have attempted to discriminate between similar species of Culicoides using multivariate analysis of continuous traits (exceptions are Pag`es et al., 2009; Augot et al., 2010). In the last 20 years a morphometric approach, based on the study of positional changes in two- or three-dimensional landmarks, has been developed, giving rise to what is termed geometric morphometrics. Practical and methodological advantages of geometric over traditional methods have revolutionized the field (Rohlf & Marcus, 1993), producing a renaissance of this discipline. Insect wings are one of the most commonly employed morphological structures in geometric morphometric studies owing to taxonomic importance and suitability for two-dimensional analyses by means of digital images. To our knowledge, however, only one work has employed this technique to analyse inter-specific wing form variation in the genus Culicoides (Mu˜noz-Mu˜noz et al., 2011). In that work significant differences between cryptic species of the subgenus Culicoides were detected, highlighting the potential of this approach to study morphological changes in other species complexes of the genus (Mu˜noz-Mu˜noz et al., 2011). In addition to species identification, the study of the phylogenetic relationships within the Obsoletus group has also been considered a relevant issue. Several previous works have attempted to reconstruct the phylogeny of the Obsoletus group employing molecular markers such as the cytochrome oxidase subunit I (COI), and internal transcribed spacer regions I (ITSI) and II (ITSII) (Meiswinkel et al., 2004; Gomulski et al., 2005; Mathieu et al., 2007; Nolan et al., 2007; Schwenkenbecher et al., 2009). Phylogenetic relationships among species of the Obsoletus group obtained from these studies are far from conclusive and conflicting results have been produced (reviewed by Schwenkenbecher et al., 2009). In the present work, wing form of the main species of the Obsoletus group and C. imicola is analysed for the first time from the perspective of geometric morphometrics with the aim of detecting differentiating features and to interpret shape changes in comparison to phylogenetic trees obtained from DNA sequence analysis. The mapping of morphological variation onto independent cladograms avoids many of the controversies inherent in the construction of phylogenies by means of morphometric variables (Klingenberg & Gidaszewski, 2010) and allows testing the congruence of shape data with phylogenetic hypotheses (Klingenberg & Gidaszewski, 2010).

Material and methods Sample composition and preparation Five different species of the subgenus Avaritia: C. obsoletus, C. scoticus, C. chiopterus, C. dewulfi and C. imicola and one species of the subgenus Culicoides: C. newsteadi N2 (sensu Pag`es et al., 2009) were analysed in the study. A total of 261 specimens (164 females and 97 males) were used, collected from 25 different sites in four European countries: Sweden, United Kingdom (U.K.), Germany and Spain (Table 1). Specimens were dissected for geometric morphometric analysis and DNA extraction. Three legs were separated for DNA extraction, whereas the remaining body parts were mounted on glass microscope slides according to the procedure previously described in Mu˜noz-Mu˜noz et al. (2011).

Species identification Specimens were initially identified under the stereomicroscope according to their pattern of wing pigmentation (Kremer, 1965). As species of the Obsoletus group cannot be differentiated solely by wing pigmentation, males were identified to species level by examination of genitalias (Kremer, 1965), whereas females were identified by species-specific PCR assays (Pag`es & Sarto i Monteys, 2005; Nolan et al., 2007). Specimens from C. newsteadi N2, previously employed in Pag`es et al. (2009) and Mu˜noz-Mu˜noz et al. (2011), were assigned to a specific category by means of species-specific PCR assays (Pag`es et al., 2009).

Morphological data acquisition Digital images of wings were obtained using a Nikon Eclipse 90i microscope equipped with a 4× Plan Fluor Nikon objective lens and a Nikon DXM 1200F camera (Tokyo, Japan). Use of left and right wings was randomized as selection of one side may bias the results in case of differential directional asymmetry between groups. Images were scaled in order to avoid variation owing to picture resolution in the analysis of wing size. A set of 13 landmarks (LMs; Fig. 1) covering the wing surface was selected and recorded in each wing using the tpsDig PC program (Rohlf, 2001).

Molecular analyses DNA was extracted from individual (adult) Culicoides of each species using a Qiagen DNeasy Blood and Tissue Kit (Qiagen, Crawley, U.K.) according to the manufacturer’s methods. Species-specific PCR assays targeting the mitochondrial COI gene were performed (Pag`es & Sarto i Monteys, 2005; Nolan et al., 2007). Sequencing and phylogenetic analyses were based on a fragment of the mtDNA COI gene. Sequences were amplified by polymerase chain reaction (PCR) using primers C1-J-1718 and C1-N-2191 as described in Pag`es et al.,

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

Culicoides wing form in the subgenus Avaritia

3

Table 1. Specific and geographical composition of the sample.

Country

Site

Coordinates

Sample (nfemales ; nmales )*

Spain

Angl`es Aramunt Bonastre Caldes de Malavella Cellera de Ter Colunga Corbera de Llobregat Garcia La Almoraima La Galera Montagut Piera Proaza Quintos de Mora St. Feliu de Buixalleu St. Iscle de Vallalta Susqueda Vilanova de la Muga Friedberg Giessen Waldsieversdorf Torestorp Tors˚as Cumbria Devon

41◦ 57 N, 42◦ 12 N, 41◦ 13 N, 41◦ 50 N, 41◦ 58 N, 43◦ 29 N, 41◦ 24 N, 41◦ 08 N, 36◦ 17 N, 40◦ 40 N, 42◦ 14 N, 41◦ 30 N, 43◦ 12 N, 39◦ 24 N, 41◦ 48 N, 41◦ 38 N, 41◦ 59 N, 42◦ 18 N, 50◦ 22 N, 50◦ 35 N, 52◦ 33 N, 57◦ 23 N, 56◦ 21 N, 55◦ 04 N, 51◦ 11 N,

Co (4; 0) Co (2; 0), Cs (1; 0), Cd (1; 0) Ci (1; 0) Cs (0; 7), Ci (7; 11) Co (2; 0) Co (2; 0), Cch (2; 3), Cd (8; 7) Ci (1; 0) Co (1; 0), Ci (1; 0) Ci (1; 3) Cn (3; 0) Co (8; 0) Co (2; 0), Ci (2; 0) Co (1; 7), Cs (1; 7), Cd (0; 6) Ci (0; 1) Cs (2; 0), Cd (0; 1) Ci (1; 0) Co (1; 0), Cs (13; 0), Cd (2; 1) Ci (1; 0), Cn (11; 5) Co (9; 0), Cs (12; 8), Cch (4; 3) Co (0; 8) Co (2; 0) Co (18; 0) Cs (15; 0) Cd (5; 0) Co (8; 5), Cs (0; 5), Cch (9; 4), Cd (0; 5)

Germany

Sweden UK

02◦ 38 E 00◦ 59 E 01◦ 25 E 02◦ 50 E 02◦ 37 E 05◦ 21 W 01◦ 55 E 00◦ 38 E 05◦ 26 W 00◦ 27 E 02◦ 36 E 01◦ 45 E 06◦ 04 W 04◦ 07 W 02◦ 34 E 02◦ 33 E 02◦ 32 E 03◦ 02 E 08◦ 55 E 08◦ 45 E 14◦ 05 E 12◦ 39 E 15◦ 49 E 02◦ 46 W 03◦ 56 W

* The

number of specimens of each species trapped in each sampling site is indicated within brackets. nfemales , number of females; nmales , number of males; Co, C. obsoletus; Cs, C. scoticus; Cd, C. dewulfi ; Ci, C. imicola; Cch, C. chiopterus; Cn, C. newsteadi N2 (sensu Pag`es et al., 2009).

Evolutionary analyses were conducted in MEGA5 (Tamura et al., 2007).

Geometric morphometrics analyses

Fig. 1. Female wing of Culicoides scoticus indicating the position of the 13 recorded landmarks (open circles). The main veins and cells are indicated. Veins: Arc, arculus; Cu, Cu1, Cu2, cubital and branches; M, M1, M2, medial and branches; R, radial; R-M, radial-median crossvein. Cells: cu, cubital; r1, r2, first and second radial.

2009. Phylogenetic analysis was inferred using the neighbourjoining method (Saitou & Nei, 1987). Evolutionary distances were computed using the Jukes–Cantor method (Jukes & Cantor, 1969). The analysis involved seven nucleotide sequences, with one representative sequence for each species analysed and an Anopheles gambiae sequence as an outgroup. Codon positions included were first + second + third + Non-coding and all ambiguous positions were removed for each sequence pair. There were a total of 477 positions in the final dataset.

In the species of the Avaritia subgenus, the wing form defined by each set of LMs was decomposed into size and shape components, which were analysed separately. Size was defined as centroid size (CS), calculated as the square root of the sum of squared distances of each landmark to the centroid of the landmarks configuration (Bookstein, 1991). An ancova with sex and species as categorical predictors and latitude as covariable was conducted to detect the effect of both factors on wing CS, controlling for potential latitudinal variation by inclusion as a covariable (Dujardin, 2008). To test for differences between pairs of species the Tukey’s HSD (Honestly Significant Difference) procedure for unequal sample size was applied (Spjotvoll & Stoline, 1973). Shape information was obtained from configurations of LMs by generalised Procrustes superimposition and projection onto tangent space (Rohlf & Slice, 1990). Allometry can be a confounding factor when shape changes are compared between species with diverging sizes and/or allometric trajectories. Thus, inter-specific allometry was assessed in the whole sample by regressing shape variables against CS. Further, a manova

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

4 F. Mu˜noz-Mu˜noz et al. with species and sex as categorical predictors and size as a covariable was conducted in the whole sample to compare allometric patterns across species and additionally to test for the effect of sex and species in wing shape variation (Gidaszewski et al., 2009). In spite of significant allometry, only the results of analyses of total shape variation were reported because: (a) a small percentage of shape changes was explained by variation in size, (b) the dependence of shape onto size differed between species and (c) the use of the residuals of the within-group regression of shape onto size provided very similar results but lower discrimination power than raw shape variables (data not shown). As sexual dimorphism was detected shape changes of the wing were studied separately in males and females. Clinal geographical variation in wing shape was tested by a manova with species and latitude as categorical and continuous predictor variables, respectively. The classification of species and the relationships among them according to wing shape were then assessed by discriminant function and canonical variate analyses. Mahalanobis distances between the centroids of species were obtained from canonical variate analyses and its P -values calculated from permutation tests (10 000 permutation rounds). With the exception of anovas and discriminant analyses, all morphometric analyses were conducted using MorphoJ (Klingenberg, 2011).

Phylogenetic signal and homoplasy analyses To test the congruence between morphometric and phylogenetic data we followed the methods described by Klingenberg and Gidaszewski (2010) and implemented in MorphoJ (Klingenberg, 2011). A principal components analysis of the covariance matrix of shape tangent coordinates was computed from the six species means. As unequal sample sizes of males and females were used, species means were conducted by averaging the mean shape of both sexes. The data points were plotted in a shape space represented by the two first principal components, which represents an optimal representation in two dimensions accounting for the maximal amount of the variation among species (Klingenberg & Gidaszewski, 2010). The congruence of morphometric data was tested with six alternative cladogram structures representing the main phylogenetic hypothesis of the Obsoletus group (Fig. 2) reconstructed on the basis of published molecular trees (Meiswinkel et al., 2004; Gomulski et al., 2005; Mathieu et al., 2007; Nolan et al., 2007; Schwenkenbecher et al., 2009). Additionally, congruence between morphometric data and the molecular phylogeny reconstructed with our own data was also tested. The first pair of trees, A and B, represent the classic hypothesis, in which C. obsoletus, C. chiopterus, C. scoticus and C. dewulfi form a monophyletic group (the so-called Obsoletus group), the sister taxa of which is C. imicola, which together form the monophyletic subgenus Avaritia. Cladograms based on ITSII support this hypothesis (Meiswinkel et al., 2004; Gomulski et al., 2005). Trees A and B only differ in which species is the sister taxon of C. chiopterus (C. obsoletus or C. scoticus). While studies based on ITSII cluster C. chiopterus

Fig. 2. Alternative phylogenies proposed for the subgenus Avaritia.

and C. scoticus, cladograms based on ITSI support the clade formed by C. obsoletus and C. chiopterus (Mathieu et al., 2007). In the second pair of trees (C and D), the subgenus Avaritia is monophyletic, but C. dewulfi is considered the sister species of C. imicola rather than being the sister species of the clade formed by the remaining species of the Obsoletus complex. Some phylogenetic analyses on molecular data support the clustering of C. dewulfi with C. imicola (Mathieu et al., 2007; Schwenkenbecher et al., 2009). Both trees, B and C, only differ in the sister taxon of C. chiopterus. Finally, trees E and F do not support the monophyly of the Obsoletus complex nor of the subgenus Avaritia. While in

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

Culicoides wing form in the subgenus Avaritia tree E C. newsteadi N2 (as a representative of the subgenus Culicoides) is the sister taxon of the clade constituted by C. obsoletus, C. chiopterus and C. scoticus, in tree F C. newsteadi N2 is the sister taxon of a clade including C. dewulfi and C. imicola. Tree E has been obtained by maximum likelihood analysis on concatenated sequences of COI, ITSI and ITSII employing Aedes aegypti as an outgroup, whereas tree F is a modification of the cladogram obtained by the same analysis without an outgroup species (Schwenkenbecher et al., 2009). Tree G corresponds to the molecular phylogeny obtained with our own data (Fig. 3). The phylogenetic trees were projected into the morphometric space using the criterion of squared-change parsimony to reconstruct the internal nodes from the shape averages of terminal taxa. The phylogenetic signal of morphometric data for each of the trees was tested by randomly permuting the shape means among the species (Klingenberg & Gidaszewski, 2010). P -values were obtained after 10 000 permutations. Homoplasy in morphometric data was quantified by computing shape consistency and shape retention indices (Klingenberg & Gidaszewski, 2010). To obtain both indices it was necessary to compute the minimal and maximal amounts of change for any unrooted tree topology adjusted with squared-change parsimony in the given morphometric data. To find the shortest unrooted bifurcating tree connecting data in the morphometric space defined by the two principal components, the FindSteinerTree programme was used (http://www.flywings.org.uk/FindSteinerTree/; Klingenberg & Gidaszewski, 2010). The maximum change was represented by the star phylogeny and was calculated as the sum of the species means from the overall mean shape in squared Procrustes distances (Klingenberg & Gidaszewski, 2010).

Results Molecular analyses Sequences were obtained at one loci (part of mtCOI gene) for all of the five Culicoides species analysed and the two outgroup species selected (C. newsteadi N2 and A. gambiae). GenBank accession numbers for them are: KF419407 (C. scoticus), KF419408 (C. obsoletus), KF419410 (C. chiopterus), KF419406 (C. dewulfi ), KF419409 (C. imicola), KF419411 (C. newsteadi N2) and DQ792577 (A. gambiae). Figure 3 shows the optimal phylogenetic tree with the sum of branch length (0.693). In the phylogenetic tree, C. obsoletus and C. scoticus are joined in a single clade, connected to C. chiopterus. There is a second clade, connected to the former group of three species, including C. dewulfi and C. newsteadi N2. Finally, C. imicola appears in a single branch connected to the other to clades, followed by A. gambiae.

Wing size analyses The ancova performed to detect the effect of species, sex and geographic clinal variation on CS highlighted that the only factor having a significant effect on wing size was the

5

Fig. 3. Neighbour-Joining tree based on COI sequences of six species of Culicoides with Anopheles gambiae as outgroup. The bootstrap values are shown next to the branches.

Fig. 4. Mean values of centroid size (CS) of males (squares) and females (triangles) of the six analysed species of the subgenus Avaritia. The vertical bars indicate 95% confidence intervals.

species (F = 29.2; d.f. = 4; P < 0.001). A significant interaction between species and sex (F = 7.6; d.f. = 4; P < 0.001), however, indicated that an inter-specific divergence of sexual size dimorphism existed. Significant differences in CS values were detected in all paired comparisons between species (P < 0.01), except in the pairs C. obsoletus - C. dewulfi and C. chiopterus - C. imicola (Fig. 4).

Wing shape analyses A significant dependence of wing shape variation with size was detected in the whole sample (P < 0.001), but only 6.4% of shape variation was explained by this parameter. The manova results confirmed the presence of allometric shape changes (significant effect of CS; Table 2) and highlighted that the dependence differed between species (significant interaction between species and CS; Table 2). Within-group pooled regression assumes equal allometric trajectories, and as a consequence the use of the residuals from pooled data confounds within- and between-group variation (McCoy et al., 2006). Considering the unequal allometric patterns among the species of the subgenus Avaritia, the small percentage of shape variation explained by size, and that the multivariate analyses using the residuals produced results that were very similar to those of complete variation, only the results of analyses of the total shape variation were reported. Sexual

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

6 F. Mu˜noz-Mu˜noz et al. Table 2. Results of the manova testing for species, sex and size (CS) effects on wing shape.

Effect

Wilk’s lambda F

d.f.effect

d.f.error

P

Species Sex CS Species × sex Species × CS Sex × CS Species × sex × CS

0.587 0.792 0.691 0.562 0.583 0.889 0.567

88 22 22 88 88 22 88

797.3 201.0 201.0 797.3 797.3 201.0 797.3

0.037 0.001 0.000 0.009 0.030 0.303 0.012

1.31 2.40 4.08 1.42 1.32 1.14 1.40

shape dimorphism and wing shape differences between species of the subgenus Avaritia were also detected (Table 2). The significant interaction between both factors indicated that intersexual differences in wing shape were species-specific. A geographical clinal variation on wing shape was observed both in males and females; however, the significant interaction between species and latitude indicated that geographical clinal variation was again species-specific (Table 3). The discriminant analyses allowed classification into species of 91.3% of the females (λwilks = 0.030; F = 8.02; d.f. = 88, 492; P < 0.001) and 92.4% of the males (λwilks = 0.009; F = 6.73; d.f. = 88, 263; P < 0.001), with percentages ranging from 80.0% to 100% depending on the species and the sex (Table 4). The first two canonical variables (CV) accounted for a similar percentage of total variation in the two sexes: 85% in females and 88% in males. CV1 explained 62% of the variance in females and 79% in males, and CV2 23% and 9%. In both sexes CV1 clearly separated C. imicola from the remaining species (Fig. 5A, B). The species of the Obsoletus group were mainly separated by CV2 in the case of females and by a combination of CV1 and CV2 in the case of males. Shape changes associated with CV1 were concordant in males and females (Fig. 5C, D). In both sexes variation in CV1 primarily involved an elongation of the distal part of the wing and a contraction of the proximal part, affecting the cubital and the radial cells (particularly in female specimens). Divergence in the relative position of some of the species of the Obsoletus group when comparing the scatterplots of males and females indicated that the magnitude of shape changes associated with CV1 among species varied between sexes. The second canonical variable represented more complex changes of wing Table 3. Results of the manovas testing for species and latitude effects on wing shape of females and males.

Sex

Effect

Females Species Latitude Species × latitude Males Species Latitude Species × latitude

Wilk’s lambda F

d.f.effect d.f.error P

0.340 0.723 0.337 0.194 0.499 0.197

88 22 88 88 22 88

1.626 1.829 1.703 1.422 2.976 1.403

473.0 119.0 473.0 243.7 61.0 243.7

0.001 0.021 0.001 0.019 0.000 0.023

Table 4. Mahalanobis distances between male (below diagonal) and female (above diagonal) species centroids (all significant at P < 0.001 level), and percentages of specimens classified in their predicted specific category in the discriminant analyses.

Co Cs Cch Cd Ci

Co

Cs

Cch

Cd

Ci

%females

%males

— 2.8 4.4 4.3 9.2

2.4 — 4.1 4.6 9.5

2.8 3.9 — 3.7 8.1

3.6 4.1 3.0 — 6.4

6.9 6.7 8.2 6.8 —

95.0 84.1 86.7 93.8 100.0

80.0 92.6 100.0 95.0 100.0

%females , percentage of females classified in their predicted specific category; %males , percentage of males classified in their predicted specific category; Co, C. obsoletus; Cs, C. scoticus; Cch, C. chiopterus; Cd, C. dewulfi ; Ci, C. imicola.

shape and slightly diverged between males and females. In males, variation in CV2 was mainly associated with an anterior shift of cubital veins and median veins junctions, whereas in females this was associated with changes in the shape of the cubital cell and the length of the radial cells. Mahalanobis distances between pairs of species were always significant at a P < 0.001 level (Table 4).

Phylogenetic and homoplasy analyses The first two principal components analysed for wing shape accounted for 52.0% and 43.9% of the variation among species means, respectively (Fig. 6). The first principal component defined a trend of complex changes of the shape that involved a broadening of the wing, a shortening of the anal and radial cells, and important displacements of median and cubital veins branches from C. newsteadi to C. imicola, with the remaining species at an intermediate position. Changes of shape defined by the second principal component mainly involved a shortening of the proximal and an elongation of the distal parts of the wing from C. obsoletus to C. newsteadi (Fig. 6). The comparison of morphological data with molecular cladograms indicated that wing shape is phylogenetically structured (Fig. 7). A significant phylogenetic signal was detected in the seven comparisons, but the lengths of the trees changed with the phylogenetic hypothesis (Table 5). The shortest tree lengths were detected for cladograms C and D, whereas trees F and G were the longest. Cladograms A, B and E had intermediate tree lengths. The minimum and maximal expected tree lengths for the particular shape data were 0.00436 and 0.0079 squared Procrustes distances, respectively. The lengths of the tested unrooted trees ranged from 0.00462 to 0.00531 (Table 6). Trees D, E and F only differed in the placement of the root and consequently the unrooted tree of the three phylogenies was exactly the same (termed D, E, F), and lengths and homoplasy indexes were identical. The lengths of the unrooted trees C and ‘D, E, F’ were close to the minimum tree length, and therefore consistency and retention indices were high, indicating a low degree of homoplasy. The lengths of unrooted trees A, B, and

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

Culicoides wing form in the subgenus Avaritia

7

Fig. 6. Tree obtained with COI sequences (own data) superimposed onto the plot of the two first principal components of the covariance matrix among species means. The percentage of variation explained by each principal component is within parentheses. The tips of the terminal branches are at the locations of species means. The outline diagrams show the shape change from the reconstructed common ancestor (grey outlines and empty circles) to the mean shape for the respective species (black outlines and solid dots) with a scale factor of 0.1 units of Procrustes distance in the positive direction.

G were considerably higher than the others and had a higher degree of homoplasy, particularly when estimated by the retention index. This index was lower than the consistency index in all the tested topologies, but especially in trees A and B. The projection of the minimal-length tree in the plot of the first two principal components (Fig. 8) demonstrated that this tree is very similar to phylogenetic hypothesis C (Fig. 7C). In fact, the topology of unrooted tree C (graph not shown) is almost identical to the minimal-length tree, although its length is higher.

Discussion

Fig. 5. Plot of the two first canonical roots (CV1 and CV2) obtained in the canonical variate analysis of females (A) and males (B). Shape changes and percentages of variance associated to each canonical root in females and males are shown respectively in (C, D). Grey outlines and open circles represent the consensus configuration of each sex, and black outlines and circles represent a shape change of 10 units of Mahalanobis distance in the positive direction along the respective CV axis.

A significant correlation between variation in wing shape and genetic distance has been previously detected in Culicoides (Mu˜noz-Mu˜noz et al., 2011). In the present study, the application of a new geometric morphometric approach has confirmed the existence of a phylogenetic signal and revealed that differences in wing shape within the subgenus Avaritia are consistent with the theory that the subgenus Avaritia is monophyletic and that C. dewulfi should be placed outside the Obsoletus group. In addition, separation of species by this method was relatively consistent, achieving positive identification in female specimens of 84.1% (C. scoticus) to 100% (C. imicola), indicating that this technique could potentially play a role in surveillance of vector species; if not as a routine tool, at least as an economic alternative to molecular techniques when precise specific identification is needed, such in pathogen detection studies (Mu˜noz-Mu˜noz et al., 2011). A key advantage of the present study lay in the use of Culicoides from a relatively wide geographical range.

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

8 F. Mu˜noz-Mu˜noz et al.

Fig. 7. Projection of alternative tree topologies onto the plot of the two first principal components of the covariance matrix among species means. The trees superimposed in (A–F) correspond to phylogenetic hypotheses presented in Fig. 2.

This allowed the potential impact of latitudinal variation on wing form, previously documented in morphometric studies (Dujardin, 2008; Mu˜noz-Mu˜noz et al., 2011; Prudhomme et al., 2012), to be evaluated. A significant effect of latitude on wing CS was not detected. While this fact does not

in itself preclude an environmental effect, it lessens the probability that observed differences in size between species are masked by an underlying environmental factor. In fact, the ordination of the species according to CS was concordant with data on wing lengths from previous taxonomic works, with

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

Culicoides wing form in the subgenus Avaritia Table 5. Tree lengths computed with unweighted squared-change parsimony and P values for the permutation test of phylogenetic signal.

Phylogeny

Tree length

P

A B C D E F G

0.00388 0.00393 0.00370 0.00379 0.00401 0.00438 0.00429

0.007 0.014 0.005 0.012 0.014 0.022 0.0428

Note: The phylogenies A–F correspond to those depicted in Fig. 2 and phylogeny G to that in Fig. 3. Rooted trees were used for these analyses. Table 6. Shape consistency and retention indices.

Phylogeny

Tree length

Consistency index

Retention index

A B C D, E, F G

0.00524 0.00531 0.00462 0.00475 0.00511

0.832 0.822 0.943 0.918 0.854

0.748 0.730 0.925 0.889 0.787

Note: The phylogenies A–F correspond to those depicted in Fig. 2 and phylogeny G that in Fig. 3. The tree lengths were computed using unweighted squared-changed parsimony and for unrooted trees.

Fig. 8. Minimal-length tree for wing shape of the six studied species.

C. scoticus and C. chiopterus having the biggest and the smallest wings, respectively, and the remaining species having wings of intermediate size (Kremer, 1965). Interestingly, a latitudinal variation in wing shape was detected, even when using size-corrected shape variables (not shown), the magnitude of which was dependent upon species. These allometry-free shape differences of environmental origin may arise between conspecific populations from distant ecogeographical areas and have been noted previously (Dujardin,

9

2008). These results indicate that changes of wing shape between distant conspecific populations of Culicoides are not a result of size variation and suggest that genetic differences may arise as species-specific adaptation to particular environments. Studies specifically addressing shape variation of conspecific populations throughout latitudinal clines are needed, however, to assess in detail the impact of latitude and underlying environmental factors, in addition to modelling the form of the wing in species of Culicoides with wide geographical distributions, such as C. obsoletus, C. scoticus or C. imicola. The impact of sexual dimorphism on the relative size of the wing varied significantly according to species. Males of C. scoticus and C. obsoletus were found to have larger wings than females, whereas similar sizes of wings were found in C. chiopterus and females of C. dewulfi and C. imicola possessed larger wings than males (Fig. 4). Sexual size dimorphism (SSD) is considered a key evolutionary feature reflecting different selection pressures acting on males and females (Lande, 1980). Given that the direction of intra-specific SSD has proved to be almost invariant (Teder & Tamaru, 2005), the divergences between the three groups of species suggest that different selection pressures could be acting in each of them. While implicit in species descriptions (e.g. Kremer, 1965), no other study to date appears to have explicitly analysed intersexual shape differences within the genus Culicoides. The present study confirms that wings of male and female Culicoides have different shapes and highlight that in the subgenus Avaritia some features of shape are common in the individuals of the same sex across all the species. In general, male Culicoides wings were elongated and narrow in comparison with the more rounded and broad wings of females (Fig. 5). Very similar results have been reported in Chironomus imicola, where sexual shape dimorphism (SShD) has been associated with different flight behaviours in males and females (McLachlan, 1986). Flight behaviours described in some species of Culicoides (e.g. Blackwell et al., 1992) are similar to those described in C. imicola (McLachlan, 1986), supporting the hypothesis that SSD and SShD could be related to different flight requirements. In addition to their obvious function in locomotion, the wings may play an important role in courtship and mate recognition in some species of Culicoides (Yuval, 2006). Accordingly, interspecific differences in sexual dimorphism of the wing could be the result of diverse selection pressures in these particular traits acting with different intensities in each sex within each species. Across all species changes in size were responsible for just 6.4% of shape variation, a value that is very close to the 5.7% observed previously in the species of the subgenus Culicoides (Mu˜noz-Mu˜noz et al., 2011). A low association between size and shape variation in Drosophila wings has been considered to be consistent with few shared quantitative trait loci determining phenotype (Zimmerman et al., 2000). In the present study, allometric patterns diverged between species, indicating that the strength of the dependence of wing shape on size is not constant in closely related species of the subgenus Avaritia.

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

10

F. Mu˜noz-Mu˜noz et al.

The shape of Culicoides wings diverged significantly between species, both in males and females. A high percentage of specimens could be assigned correctly to species according to wing shape (more than 90%). Culicoides imicola possessed the most divergent wing shape within the subgenus Avaritia, with no overlap in the shape space with the other species analysed. The species of the Obsoletus group possessed more similar wing shapes and some inter-specific overlap existed, with the exception of C. chiopterus males. The high percentage of discrimination detected in the species of the Obsoletus group, however, supported the idea that the multivariate study of wing shape is an effective tool to distinguish cryptic species of Culicoides (Mu˜noz-Mu˜noz et al., 2011). Although the power of discrimination between species of the Obsoletus group obtained by multivariate analysis of body measurements has proven somewhat higher in specimens trapped in France (see Augot et al., 2010), the identification of specimens of cryptic species by means of wing shape can be useful, particularly if the remaining body of the Culicoides is required for further analyses (e.g. for pathogen detection: Mu˜nozMu˜noz et al., 2011). The geometric morphometric data were notable for high consistency and retention indices obtained during the seven phylogenetic comparisons, particularly in the case of trees C, D, E and F. This indicated a low degree of homoplasy in wing shape of the Culicoides examined and less than previously detected in Drosophila species (Klingenberg & Gidaszewski, 2010). When trees D, E and F were rooted, the lengths of the trees E and F were somewhat higher than that of D. These results imply that wing shape is consistent with monophyly of the subgenus Avaritia. In previous studies based on molecular markers, phylogenetic analyses of ITSI sequences of several species of Culicoides from France supported the existence of a monophyletic clade formed by some of the species of the subgenus Avaritia (C. obsoletus, C. scoticus and C. imicola; Perrin et al., 2006), however, the simultaneous analysis of COI, ITSI and ITSII sequences of a more restricted group of species of Culicoides rejected such a clade (Schwenkenbecher et al., 2009). In addition, the neighbourjoining tree calculated with our own data clustered together C. newsteadi N2 and C. imicola, rejecting the monophyly of the subgenus Avaritia. Resolving this issue will require the use of a wider range of species sourced from disparate locations as part of an ongoing reassessment of subgenus classification in Culicoides. Among the seven tested phylogenetic hypotheses, C was the most concordant with wing shape. This phylogeny not only had the shortest tree length and the highest consistency and retention indices, but the projection of the unrooted version of the tree was almost identical to the minimal treelength topology obtained by the Steiner method. Thus, shape variation of the wing is also consistent with the exclusion of C. dewulfi from the Obsoletus group, previously proposed after the simultaneous phylogenetic analysis of multiple molecular markers (Schwenkenbecher et al., 2009). A third phylogenetic inference lies in the relationship between C. obsoletus, C. scoticus and C. chiopterus. Most phylogenetic studies concur in that these three species form a clade with C. montanus in regions where this relatively

rare species is present (Meiswinkel et al., 2004; Gomulski et al., 2005; Mathieu et al., 2007). Evolutionary relationships among these species remain unclear, however, and three possible assemblages have been obtained in different studies. Molecular phylogenies obtained with COI (Nolan et al., 2007) and concatenated sequences of COI, ITSI and ITSII (Schwenkenbecher et al., 2009) grouped C. scoticus with C. chiopterus; a phylogenetic tree using ITSI sequences clustered C. obsoletus (and C. montanus) with C. chiopterus (Mathieu et al., 2007); and the phylogenetic analyses of COI sequences from the present study grouped C. obsoletus and C. scoticus as sister species. While C. obsoletus and C. scoticus have the most similar wing shapes among the species of the subgenus, both in males and females, the tree of minimal length clustered together C. obsoletus and C. chiopterus, agreeing with the ITSI phylogeny. Moreover, cladograms with this grouping (phylogenies A and C) had lower tree lengths than those that clustered C. scoticus and C. chiopterus (phylogenies B and D). The present study has demonstrated that geometric morphometrics provides a useful tool for separation of cryptic Culicoides of importance as vectors. While the results of wing analysis largely were concordant with traditional taxonomic groupings, established mainly by qualitative traits regarding wing markings and genitalia, important discrepancies also appeared. The most important of these was the exclusion of C. dewulfi from the Obsoletus group suggested by means of wing shape, a conclusion not originally drawn from morphological studies. As a whole, the study highlights a need for a consistent and systematic phylogeny of the subgenus, which integrates morphological traits and several molecular markers (including additional types of data, such as karyotypes, if possible). Finally, the results underline that the morphological discrimination between cryptic species is not only a consideration of how variable they are but of the precision of the devices and the techniques that we use to distinguish them. In this sense, geometric morphometrics has the potential to be a powerful tool to taxonomists and other researchers in discriminating among cryptic species of Culicoides.

Acknowledgements The authors would like to thank the owners of farms for their permission to place the traps there and people involved in sample collection. We also thank Departament d’Agricultura, Alimentaci´o i Acci´o Rural (DAR) for financing the BT Entomological Surveillance Network and for logistic support. This study was partly financed by the Instituto Nacional de Investigaciones Agrarias (project number FAU2008-00019C03). S. Carpenter was funded by the Defra Culicoides Reference Laboratory, U.K. References Augot, D., Sauvage, F., Jouet, D. et al. (2010) Discrimination of Culicoides obsoletus and Culicoides scoticus, potential bluetongue vectors, by morphometrical and mitochondrial cytochrome oxidase subunit I analysis. Infection, Genetics and Evolution, 10, 629–637.

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042

Culicoides wing form in the subgenus Avaritia Blackwell, A., Mordue, L., Young, M.R. & Mordue, W. (1992) The swarming behaviour of the Scottish biting midge, Culicoides impunctatus (Diptera: Ceratopogonidae). Ecological Entomology, 17, 319–325. Bookstein, F.L. (1991) Morphometric Tools for Landmark Data. Cambridge University Press, Cambridge, U.K. Borkent, A. & Wirth, W.W. (1997) World species of biting midges (Diptera: Ceratopogonidae). Bulletin of the American Museum of Natural History, 233, 1–257 (updated 2011). Cˆetre-Sossah, C., Baldet, T., Del´ecolle, J.C., Mathieu, B., Perrin, A., Grillet, C. & Albina, E. (2004) Molecular detection of Culicoides spp. and Culicoides imicola, the principal vector of bluetongue (BT) and African horse sickness (AHS) in Africa and Europe. Veterinary Research, 35, 325–337. Dujardin, J.P. (2008) Morphometrics applied to medical entomology. Infection, Genetics and Evolution, 8, 875–890. Gidaszewski, N.A., Baylac, M. & Klingenberg, C.P. (2009) Evolution of sexual dimorphism of wing shape in the Drosophila melanogaster subgroup. BMC Evolutionary Biology, 9, 110–120. Gomulski, L.M., Meiswinkel, R., Del´ecolle, J.C., Goffredo, M. & Gasperi, G. (2005) Phylogenetic relationships of the subgenus Avaritia Fox, 1955 including Culicoides obsoletus (Diptera, Ceratopogonidae) in Italy based on internal transcribed spacer 2 ribosomal DNA sequences. Systematic Entomology, 30, 619–631. Hoffmann, B., Scheuch, M., Hoper, D. et al. (2012) Novel orthobunyavirus in cattle, Europe, 2011. Emerging Infectious Diseases, 18, 469–472. Jukes, T.H. & Cantor, C.R. (1969) Evolution of protein molecules. Mammalian Protein Metabolism (ed. by H.N. Munro), pp. 21–132. Academic Press, New York, NY. Klingenberg, C.P. (2011) MorphoJ: an integrated software package for geometric morphometrics. Molecular Ecology Research, 11, 353–357. Klingenberg, C.P. & Gidaszewski, N.A. (2010) Testing and quantifying phylogenetic signals and homoplasy in morphometric data. Systematic Biology, 59, 245–261. Kremer, M. (1965) Contribution a` l’´etude du genre Culicoides Latreille particuli`erement en France, Encyclop´edie Entomologique, S´erie A, ´ Vol. 39 . Editions Paul Lechevalier, Paris. Lande, R. (1980) Sexual dimorphism, sexual selection and adaptation in polygenic characters. Evolution, 34, 292–305. Mathieu, B., Perrin, A., Baldet, T., Del´ecolle, J.C., Albina, E. & C´etreSossah, C. (2007) Molecular identification of Western European species of obsoletus complex (Diptera: Ceratopogonidae) by an internal transcribed spacer-1 rDNA multiplex polymerase chain reaction assay. Journal of Medical Entomology, 44, 1019–1025. McCoy, M.W., Bolker, B.M., Osenberg, C.W., Miner, B.G. & Vonesh, J.R. (2006) Size correction: comparing morphological traits among populations and environments. Oecologia, 148, 547–554. McLachlan, A.J. (1986) Sexual dimorphism in midges: strategies for flight in the rain-pool dweller Chironomus imicola (Diptera, Chironomidae). Journal of Animal Ecology, 55, 261–267. Meiswinkel, R., Gomulski, L.M., Del´ecolle, J.C., Goffredo, M. & Gasperi, G. (2004) The taxonomy of Culicoides vector complexes – unfinished business. Veterinaria Italiana, 40, 151–159. Mellor, P.S., Boorman, J.P.T. & Baylis, M. (2000) Culicoides biting midges: their role as arbovirus vectors. Annual Review of Entomology, 45, 307–340.

11

Mu˜noz-Mu˜noz, F., Talavera, S. & Pag`es, N. (2011) Geometric morphometrics of the wing in the subgenus Culicoides (Diptera: Ceratopogonidae): from practical implications to evolutionary interpretations. Journal of Medical Entomology, 48, 129–139. Nielsen, S.A. & Kristensen, M. (2011) Morphological and molecular identification of species of the Obsoletus group (Diptera: Ceratopogonidae) in Scandinavia. Parasitology Research, 109, 1131–1141. Nolan, D.V., Carpenter, S., Barber, J., Mellor, P.S., Dallas, J.F., Mordue, L. & Piertney, S.B. (2007) Rapid diagnostic PCR assays for members of the Culicoides obsoletus and Culicoides pulicaris species complexes, implicated vectors of bluetongue virus in Europe. Veterinary Microbiology, 124, 82–94. Pag`es, N. & Sarto i Monteys, V. (2005) Differentiation of Culicoides obsoletus and Culicoides scoticus (Diptera: Ceratopogonidae) based on mitochondrial cytochrome oxidase subunit I. Journal of Medical Entomology, 42, 1026–1034. Pag`es, N., Mu˜noz-Mu˜noz, F., Talavera, S., Sarto, V., Lorca, C. & N´un˜ ez, J.I. (2009) Identification of cryptic species of Culicoides (Diptera: Ceratopogonidae) in the subgenus Culicoides and development of species-specific PCR assays based on barcode regions. Veterinary Parasitology, 165, 298–310. Perrin, A., Cetre-Sossah, C., Mathieu, B., Baldet, T., Del´ecolle, J.C. & Albina, E. (2006) Phylogenetic analysis of Culicoides species from France based on nuclear ITS1-rDNA sequences. Medical and Veterinary Entomology, 20, 219–228. Rohlf, F.J. (2001) tpsDig, Version 1.31 . State University, New York, NY. Rohlf, F.J. & Marcus, L.F. (1993) A revolution in morphometrics. Trends in Ecology and Evolution, 8, 129–132. Rohlf, F.J. & Slice, D. (1990) Extensions of the procrustes method for the optimal superimposition of landmarks. Systematic Zoology, 39, 40–59. Saegerman, C., Berkvens, D. & Mellor, P.S. (2008) Bluetongue epidemiology in the European Union. Emerging Infectious Diseases, 14, 539–544. Saitou, N. & Nei, M. (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular Biology and Evolution, 4, 406–425. Schwenkenbecher, J.M., Mordue, A.J. & Piertney, S.B. (2009) Phylogenetic analysis indicates that Culicoides dewulfi should not be considered part of the Culicoides obsoletus complex. Bulletin of Entomological Research, 99, 371–375. Spjotvoll, E. & Stoline, M.R. (1973) Extension of t-method of multiple comparison to include cases with unequal sample sizes. Journal of the American Statistical Association, 68, 975–978. Tamura, K., Dudley, J., Nei, M. & Kumar, S. (2007) MEGA4: molecular evolutionary genetics analysis (MEGA) software version 4.0. Molecular Biology and Evolution, 24, 1596–1599. Teder, T. & Tamaru, T. (2005) Sexual size dimorphism within species increases with body size in insects. Oikos, 108, 321–334. Yuval, B. (2006) Mating systems of blood-feeding flies. Annual Review of Entomology, 51, 413–440. Zimmerman, E., Palsson, A. & Gibson, G. (2000) Quantitative trait loci affecting components of wing shape in Drosophila melanogaster. Genetics, 155, 671–683. Accepted 16 September 2013

© 2014 The Royal Entomological Society, Medical and Veterinary Entomology, doi: 10.1111/mve.12042