Sample Collection

Samples of Wilhelmia species were collected from 2014 to 2017 at 21 sites in the Balkan Peninsula (Table 1, Fig. 1). All samples (larvae and pupae) were fixed in 96% ethanol. The individuals were identified morphologically to the species level when possible or to the lowest possible species group level using different identification keys (Rubtsov 1956, Knoz 1965, Rivosecchi 1978, Crosskey 2002, Yankovsky 2002, Jedlička et al. 2004, Lechthaler and Car 2005).

Isolation of DNA, PCR Amplification, and Sequencing

Genomic DNA was extracted from each individual simuliid using the KAPA2G Express Extract Kit following the manufacturer's instructions. The quality of the DNA was tested on a 1% agarose gel. For 47 individuals of six morphologically identified species—S. balcanicum (18), S. turgaicum (2), S. lineatum (3), S. pseudequinum (9), S. equinum (11) and S. paraequinum (4)—the barcoding region of the mitochondrial COI gene was amplified using primers: LCO1490 (5-GGTCAACAAATCATAAAGATATTGG-3) and HCO2198 (5-TAAACTTCAGGCTGACCAAAAAATCA-3) (Folmer et al. 1994). mtDNA amplification was performed twice in a volume of 25 μl. The reaction mixture contained 1 μl of extracted DNA, 16.9 μl of dH₂O, 0.5 μl dNTPs, 0.5 μl GoTaq buffer, 0.7 μl of both primers, and 0.2 μl of GoTaq polymerase. PCR cycles were performed using the following thermal profiles: initiation of denaturation at 95°C for 2 min, followed by 35 cycles: 1 min of denaturation at 94°C, 1 min of primer annealing at 50°C and 1 min extension at 72°C, and the final extension step at 72°C for 5 min. PCR products were visualized on 1% agarose gels stained with ethidium bromide. DNA sequencing was performed at the Center for Human Molecular Genetics at the Faculty of Biology, University of Belgrade. All sequences were checked and arranged using the ABI Sequence Scanner Software v. 2.0 (Applied Biosystems). DNA sequences were archived at GenBank, under the accession numbers shown in Table 1.

Genetic and Phylogenetic Analyses

Since simuliid larvae can be difficult to identify, we used the BLASTn algorithm to search for similar sequences in the GenBank database that contains unidentified or misidentified simuliids. Two obtained sequences (GenBank MF458827 and MF458826) were similar to our S. pseudequinum sequences, so we included them in the analyses. Furthermore, sequences originating from species of the Simulium (Simulium) jenningsi group were more than 90% identical to the Wilhelmia subgenus in BLAST algorithms, so we used one of the species of this group, Simulium notiale Stone and Snoddy, 1969, as an outgroup for the subgenus Wilhelmia. In total, 226 sequences were analyzed: 47 sequences of Wilhelmia species collected from the Balkan Peninsula, 47 sequences of Wilhelmia downloaded from the Barcode of Life Data System (BOLD System,  http://www.boldsystems.org/; Ratnasingham and Hebert 2007), 124 sequences of Wilhelmia species from GenBank, and 8 sequences from the GenBank database were used as outgroups: 4 S. notiale, 2 Culicoides anophelis Edwards, 1922, and 2 Thaumalea testacea Ruthe, 1831. Sequences downloaded from BOLD and GenBank are listed in Supp Table S1 (online only). The map of localities from which the Wilhelmia sequences were obtained is shown in Fig. 1. All sequences were aligned in MEGA6 (Tamura et al. 2013) with the ClustalW algorithm.

Maximum likelihood (ML) and maximum parsimony (MP) phylogenetic analyses were performed using MEGA6 software (Tamura et al. 2013). We found the best-fitting models of sequence evolution in MEGA6 with the model comparison by Bayesian information criterion (BIC) and log-likelihood (lnL) and used them in the subsequent analyses. One thousand bootstrap replicates were performed to assess branch support in the resulting ML and MP trees. The best-fitting model of base substitution was also used to calculate the average genetic distances between the sequences within each clade and between clades of species by the bootstrap method (1,000 replicates) in MEGA6.

Bayesian phylogenetic analyses were performed in BEAST v2.4.2 (Bouckaert et al. 2014). For site evolution model priors, the best fitting one out of the available models within BEAST was selected according to the model selection run in MEGA6. We ran preliminary tests to examine the performance of strict versus uncorrelated log-normal relaxed clock priors. These preliminary analyses consisted of two independent runs, each for 6,000,000 generations, with sampling every 1,000 generations. We examined posterior density histograms in TRACER v1.6 (Rambaut et al. 2014), concluded that strict clock priors better model our data, and used these clock priors in constructing final gene trees.

Temporal patterns of diversification in Wilhelmia were also explored by estimating the divergence times using the Bayesian Markov chain Monte Carlo method with strict clock priors. As no fossils are available to calibrate the nodes within Wilhelmia, we used the estimates of Bertone et al. (2008) for two nodes, MRCA for Ceratopogonidae and Simuliidae+Thaumaleidae within Culicomorpha, and MRCA1 for Thaumaleidae and Simuliidae divergence. For these nodes (constrained to be monophyletic), we enforced log-normal priors with means of 226 and 130 million years ago (Ma), respectively. The important branching events in evolutionary history of Wilhelmia were dated, and the corresponding paleogeographical maps were consulted (Kazmin and Natapov 1998) in order to describe the cladogenesis and speciation events.

Table 1.

Data for Wilhelmia species collected in 2014–2017

Fig. 1.

Localities of the collected Wilhelmia specimen sequences from the Balkan Peninsula and localities of the origin for the downloaded sequences from NCBI GenBank and BOLD System.

The Bayesian phylogenetic analysis was independently run twice for 10 million generations each, with sampling every 1,000 generations. The results were visualized in Tracer v1.6 (Rambaut et al. 2014) to assess the convergence by effective sample sizes (ESS). In all cases, a burn-in of 10% was appropriate. We used LogCombiner v1.8.2. (Rambaut and Drummond 2015) to combine independent runs, with 10% of each sample discarded as burn-in and Tree annotator (Drummond et al. 2012) to generate a maximum clade credibility tree based on the mean height of clades in the posterior distribution.

Evolutionary significant units in the summarized phylogenetic tree were delimited by the single-threshold Generalized Mixed Yule Coalescent (GMYC) method (Pons et al. 2006, Fujisawa and Barraclough 2013). This method was shown to be robust for species delimitation when using single locus data (Fujisawa and Barraclough 2013). We compared the results from GMYC to morphological identification of specimens and chose to refer to the recovered units as clades rather than species.

Nucleotide diversity calculations and tests of neutrality were performed for each clade (potential species) delimited by GMYC analysis in DnaSP v6.10.01 (Rozas et al. 2017). The following parameters were determined: number of used sequences (n), number of haplotypes (h), number of segregating sites (S), haplotype diversity (Hd) with the SD, nucleotide diversity (Pi) with the SD, Tajima's D statistic (Tajima 1989), and Fu's Fs (Fu 1997). The mismatch distribution tests were performed with the program DnaSP v6.10.01 to test for changes in population size. The network of Wilhelmia haplotypes recognized in DnaSP was constructed in Network v5.0.0.1. (Librado and Rozas 2009). First, the star contraction (Forster et al. 2001) of haplotypes was conducted to reduce the number of nodes in the network. Afterwards, the median-joining algorithm (Bandelt et al. 1999) was used for network calculation.

Equinum Branch

The diagnostic karyological features of the equinum branch, according to previous studies (Weber and Grunewald 1989, Chubareva et al. 2007; Huang et al. 2012, Adler et al. 2015), would be the following: chromosomes are disconnected (not creating a chromocenter), an expanded region is present in chromosome I, and the bulge marker and ring of Balbiani are in the more distal, subterminal region of IIS.

Fig. 2.

Bayesian phylogenetic tree based on the COI gene of Wilhelmia specimens, with S. notiale, Culicoides anophelis, and Thaumalea testacea as the outgroups. The numbers above/below the branches represent posterior BA probabilities followed by ML and MP >50% bootstrap support. The sequences are given as GenBank or BOLD accession numbers. Sequences in bold type were obtained in this study. The bar shows the number of substitutions.

Within the equinum branch, the samples were morphologically identified as belonging to two widespread species, S. equinum and S. pseudequinum, that were subsequently shown to be non-monophyletic. Six clades in the phylogenetic tree (Fig. 2) were confirmed by both GMYC and the haplotype network. The genetic distance between the two clades (S. pseudequinum W1 and S. pseudequinum W2) was low (2.59%), which is in line with the most recent divergence between them. The other distances among clades in this phylogenetic branch were high (8.33–17.79%), which implies a possible species level of each single clade.

The majority of S. equinum sequences formed the monophyletic clade we named S. equinum. This youngest clade in the branch is also the most widespread and haplotypically the most diverse. A small early diverging clade of ‘S. equinum’ sequences from Turkey and Finland was named S. equinum 2. The taxonomic character for S. equinum are swollen gills without a basal constriction (Rivosecchi 1978), observed in pupal samples from Finland (BOLD System). However, several subspecies of this species were described in the keys for eastern Europe (Rubtsov 1956, Yankovsky 2002), sometimes as different species (e.g., Simulium ivashentzovi Rubtsov, 1940). Samples from all the subspecies are needed in order to better define this cryptic clade recovered from Turkey and Finland and to look for any morphological traits that can discriminate it from S. equinum.

Table 3.

Nucleotide diversity calculations and tests of neutrality

Samples identified as S. pseudequinum were grouped in four diverging clades. The calculated timing of the differentiation between the clades was between 46 and 21 Ma. At that time, the regions now harboring these clades were continents and islands which were well separated by the sea that might have allowed for allopatric speciation between the clades to happen. The S. equinum clade was nested among these clades, thus making the widespread taxon S. pseudequinum paraphyletic. Although it had never been tested before, Cherairia et al. (2014) suggested that S. pseudequinum could be a species complex due to its diverse habitats and broad geographical range. The species S. pseudequinum (Séguy, 1921) was originally described from the Canary Islands that are geographically closest to the origin of the S. pseudequinum W2 samples. Therefore, the sequences from the Canary Islands population would be of utmost importance in resolving the nomenclature within the equinum branch.

In S. pseudequinum, some authors have insisted on recognizing morphologically different subspecies, such as Wilhelmia mediterranea sulfuricola Rivisecchi, 1972 and Wilhelmia mediterranea fluminicola Rivosecchi, 1972 (Rivosecchi, 1978). Furthermore, at the beginning of the 20th century, a new species was described from Serbia, Simulium brnizense Baranov, 1924, which was later synonymized with S. pseudequinum. We assumed that the clade B of S. pseudequinum might correspond to S. brnizense. However, additional morphological, cytogenetic, and genetic analyses, which would also include type specimen of S. pseudequinum from Canary Island, the Balkans and the Middle East, need to be performed in order to resolve these issues.

Lineatum Branch

The diagnostic karyological features of our lineatum branch would be as follows (Petrova et al. 2003, Chubareva et al. 2007, Huang et al. 2012, Adler et al. 2015): chromocenter present (variable extent observed in different species), expanded region of chromosome I absent, and the bulge marker and ring of Balbiani not subterminal in IIS. Furthermore, a higher inversion polymorphism is noted in the species of the lineatum branch. In this study, three monophyletic lineages were observed within the lineatum branch. The sergenti lineage was represented by a single clade, the paraequinum lineage consisted of two clades, while the lineatum lineage was comprised of six clades. The differences between the latter groupings were previously reported, and segregation of the groupings was detected by the position of the Balbiani ring (closer to the centromere in S. paraequinum) and the size of the chromocenter (Petrova et al. 2003, Chubareva et al. 2007; Adler et al. 2015).

Within the subgenus Wilhelmia, sergenti lineage is a morphologically distinct. It is characterized by pupal gills with only four central filaments: two external filaments are long and sturdy and two dorsal filaments are short and flexible (Rubtsov 1956). The sergenti lineage morphologically includes western Palearctic species, S. sergenti and Simulium quadrifila Grenier, Faure & Laurent, 1957, but also Simulium xingyiense Chen & Zhang, 1998 of the eastern Palearctic. Although the species was not karyologically surveyed, S. xingyiense showed characteristics similar to other lineages of the lineatum branch (Huang et al. 2012). Molecular sequences included in this study represented a single clade and the species S. sergenti, and revealed a huge divergence from other clades (41 mutational steps in the haplotype network, a range of genetic distances 12.12–16.13 %).

There were two clades of S. paraequinum in all phylogenetic trees. The divergence between these clades occurred 23–11 Ma. The genetic difference between them was relatively high (5.10%), even though they were sampled from geographically neighboring localities. Two syntopical varieties (Wilhelmia paraequina paraequina Puri and Wilhelmia paraequina transcaucasica Rubzov) of S. paraequinum were first described by Rubtsov (1956) in Armenia. These varieties were distinguished by the size of the adult flies and the differences in morphology of (adult male) gonofurca. Later, two syntopical cytoforms (cytotypes A and B) of this species were described from Armenia by Petrova et al. (2003), and two mtDNA lineages were revealed in Turkey by Inci et al. (2017). Our tree shows the ancient mtDNA divergence of S. paraequinum, which could be interpreted as the existence of two taxa (possibly matching with paraequinum and transcaucasicum, per Rubtsov). As we did not have samples from Iran of the species Simulium lurestanicum (Yankovsky, 2010) that are thought to be a synonym of S. paraequinum by Khazeni et al. (2013), we could not suggest the identity nor the position of that taxon.

Within the lineatum lineage, the sibling species, S. lineatum, S. turgaicum, and S. balcanicum were present with two clades each. This is corroborated with the results of GMYC according to which the lineatum lineage harbors six species. All clades within the lineage had a low genetic distance from each other (1.78–3.30%), and the phylogenetic relationships of some clades on the trees were less conclusive (BI < 0.8). Inci et al. (2017) also found a low genetic distance between the species of the lineatum lineage (2.7–3.4 %). The divergence between these clades was shown to be recent, with the time of the earliest divergence being for lineatum ∼20–6 Ma, while the time of other divergences largely overlapped. Only within this lineage, the branches of the haplotype network highly overlapped, which further raises doubt that the clades represent different species. A similar position is obtained in the BOLD System where all six clades correspond to one cluster (BOLD: AAM4036).

Table 4.

Evolutionary divergence between clades based on the pairwise analysis of COI sequences calculated by the Tamura three-parameter method using MEGA 6 (Tamura et al. 2013)

The aforementioned species lack reliable morphological traits to be distinguished by. Some studies from Central Europe have always doubted the validity of the differences between S. lineatum and S. balcanicum, mostly when adults were compared (Crosskey and Zwick 2007, Jedlička and Seitz 2008). In some identification keys, S. turgaicum is characterized by a greater length of the anteriormost filaments of the pupal gill than S. lineatum (Rubtsov 1956); however, Adler et al. (2015) did not consider this diagnostic feature as adequate and claimed that these two species were morphologically indistinguishable. Thus, among the species of the lineatum lineage, the only obvious morphological character available is the presence of a petiolate (forked) pair of gill filaments in S. balcanicum pupae (Rubtsov 1956, Yankovsky 2002, Lechthaler and Car 2005, Jedlička and Seitz 2008, Adler et al. 2015). Our specimens that had petiolate gills and were identified as S. balcanicum were positioned in two different clades of the phylogenetic tree. One clade (S. balcanicum), which included our and downloaded sequences, was comprised of pupae and larvae that were unequivocally identified as S. balcanicum. The other clade (S. balcanicum/turgaicum 1) was comprised of sequences from Balkan individuals with clearly visible petiolate gills (therefore, S. balcanicum), but also of sequences from Turkish individuals morphologically and cytogenetically identified as S. turgaicum 1 (Inci et al. 2017). This should question petiolate gills as a reliable morphological character to distinguish between morphospecies S. balcanicum and S. turgaicum.

Adler et al. (2015) in their study recognized (chromosomally) four separate species from the lineatum grouping (S. balcanicum, S. lineatum, S. takahasii, and S. turgaicum). We could not obtain sequences from S. takahasii and will not discuss the position of this species. The chromosomal relationships of three other species given therein (S. turgaicum as sister to S. lineatum and S. balcanicum) contradict our phylogenetic scenario. The fact that S. balcanicum was shown to be very monomorphic and to bear a possibly derived and almost fixed character of IL-14 inversion (Adler et al. 2015) could be interpreted within the limits of our scenario if we assume that it shares with the clade S. turgaicum 1 the same inversion. If this turns out to be an accurate assumption, the inversion IL-14 would be a synapomorphy for clades S. balcanicum, S. balcanicum/turgaicum 1, and S. turgaicum 1. Further investigation within S. turgaicum 1 is needed to confirm this possibility. Sequences that are found scattered in S. turgaicum 1 and S. balcanicum/turgaicum 1 were obtained from the study of Inci et al. (2017), where the authors reported 100% ambiguous identification within this species. The clade S. turgaicum 2, which could then be seen as S. turgaicum s.str., is more widely distributed than thought (it is present in the Balkan peninsula as well) and should also be chromosomally investigated from the new localities to check for the reported unique states in IIL-8, IIL-11, and IIL-12 inversions (Adler et al. 2015). Two recognized clades of S. lineatum species largely overlap in distribution, as sequences from both clades were found in neighboring localities. Chromosomal differences were encountered in English specimens studied by Adler et al. (2015), where the authors observed seven individuals with linked inversions in IS and in IIL. We agree with Adler et al. (2015) that the possibility of cryptic taxa in S. lineatum require additional sampling.

Fig. 3.

Haplotype network obtained from Wilhelmia species mtCOI gene sequences using Network (Librado and Rozas 2009). Circle sizes are proportional to the haplotype frequency.

Wide distributions of recognized species groupings within the lineatum branch (west Palearctic sergenti, quadrifila + eastern Palearctic xingyiense; west/central Palearctic lineatum, balcanicum, turgaicum + eastern Palearctic takahasii) call for more sampling and a thorough phylogenetic study of central and eastern Palearctic species in order to obtain the complete scenario of Wilhelmia evolution.