Phylogeography of Simulium Subgenus Wilhelmia (Diptera: Simuliidae)-Insights From Balkan Populations.

Many morphologically similar species of the simuliid (Diptera: Simuliidae) subgenus Wilhelmia, Enderlein are difficult to distinguish. Thus, the revision of the subgenus using various morphological, cytogenetic, and genetic analyses has been attempted. Neglected until now, the Balkan Peninsula, a crossroad between Europe and Anatolia, provides insight which could resolve problematic interrelationships of the taxa within this subgenus. To uncover the status and relations within the subgenus Wilhelmia, mtDNA was extracted from 47 individuals of six morphospecies: Simulium balcanicum (Enderlein, 1924), Simulium turgaicum Rubtsov, 1940, Simulium lineatum (Meigen, 1804), Simulium pseudequinum Séguy, 1921, Simulium equinum (Linnaeus, 1758), and Simulium paraequinum Puri, 1933 from 21 sites throughout the Balkan Peninsula. Phylogenetic analysis of the Wilhelmia species using mitochondrial DNA barcoding (COI) gene showed two major branches, the lineatum branch, which includes the lineages sergenti, paraequinum, and lineatum, and the equinum branch. In the equinum branch, the mtDNA sequences formed six clades, with high genetic distances, suggesting the existence of different species. Historically, the clades of the equinum branch appeared at numerous islands, perhaps as a result of allopatric speciation. The paraequinum lineage (lineatum branch) is composed of two species. However, six clades of the lineatum lineage overlapped with intra- and interspecific genetic distances. Our results revealed that the species S. balcanicum, S. pseudequinum B, and S. equinum were omnipresent in the Balkans. The results point to not only the fair diversity of Wilhelmia species in the Balkans, but also indicate that most Wilhelmia species live in sympatry.

Black flies (Simuliidae) are cosmopolite holometabolous insects. The development of early stages (egg, larva, and pupa) is bound to running water (Currie and Adler 2008, Low et al. 2014), while adults are small flies that are considered as highly mobile, as they can cross more than 500 km (Crosskey 1990) and have even successfully colonized distant islands (Craig et al. 2001). According to the last inventory, a total of 2,351 species has been recorded (Adler and Crosskey 2018). Nevertheless, new species are identified every day, partly because of the use of new identification techniques which have disentangled their considerable morphological similarity, and partly because we have acquired a grasp of the ongoing speciation events. Speciation is a continuous process; hence contemporary population studies show current evolutionary stages. Consequently, at the time of the analysis, the speciation boundaries were not always clear (Mallet 2008, Hendry et al. 2009, Conflitti et al. 2017).

The Palearctic subgenus Wilhelmia Enderlein belongs to the Simulium Latreille genus and includes species that are considered significant and widespread pests of humans and livestock (Crosskey 1990, Werner and Adler 2005, Sarıözkan et al. 2014, Inci et al. 2017). At present, 31 extant species of Wilhelmia are described. No fossils of the subgenus are known. Wilhelmia’s areal extends from Spain and the British Isles, throughout Europe and the Middle East, parts of Kashmir, China, and western Pakistan to Japan. On the African continent, they can be found exclusively in the north, from the Canary Islands and Morocco, through Tunisia to Libya (Crosskey 1969, Adler et al. 2015, Adler and Crosskey 2018). Due to small differences in morphological characters, misidentification and incorrect description of new species has occurred quite often. Therefore, many described species have ended up synonymized (Adler and Crosskey 2018). As a result of erroneous identification, the species distribution is questionable. Moreover, inaccurate identification could result in inadequate control measures with negative socio-economic outcomes (Hernández-Triana et al. 2012).

Novel studies of the giant polytene chromosomes extracted from the larval salivary glands have revealed a high cryptic biodiversity within the Simuliidae family (Petrova et al. 2003, Adler and Crosskey 2015). Certain taxa, once considered a single widespread species, have shown to be groupings of sister species (Rothfels 1979, Adler et al. 2010, 2015). Within one such grouping in subgenus Wilhelmia, Adler et al. (2015) clarified the cytological identity of morphologically indistinguishable species and their ranges, showing the presence of four species, Simulium balcanicum (Enderlein, 1924), Simulium lineatum (Meigen, 1804), Simulium takahasii (Rubtsov, 1962), and Simulium turgaicum Rubtsov, 1940. The identification of morphologically similar species with molecular tools also proved to be equally successful. In the majority of studies, DNA barcoding using the mitochondrial cytochrome c oxidase subunit I (COI) gene is often used to reveal the cryptic diversity of black flies (Hernández-Triana et al. 2012, Pramual and Nanork 2012, Conflitti et al. 2013, Sriphirom et al. 2014, Inci et al. 2017).

In past research that relied only on morphology, the high diversity of the subgenus Wilhelmia was reported for the Balkan Peninsula. At least 10 species were described from the Balkan Peninsula, but most have been synonymized. According to the latest revised species inventory (Adler and Crosskey 2018), six morphologically distinct species of Simulium, subgenus Wilhelmia, were reported for the Balkan Peninsula: Simulium angustifurca (Rubtsov, 1956), Simulium equinum (Linnaeus, 1758), Simulium paraequinum Puri, 1933, Simulium pseudequinum Séguy, 1921, S. balcanicum, and S. lineatum.

The aims of this study were to investigate the species diversity and distribution of the subgenus Wilhelmia in the Balkans, and to determine the positions of the specimens from the Balkans within the phylogeographical frame of the subgenus.

Materials and Methods

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).

Table 1.

Data for Wilhelmia species collected in 2014–2017

Accession numbers	Species	Stage	River	Country	Latitude	Longitude	Collection date	Collector
MH549547	Simulium balcanicum/turgaicum 1	Pupa	Lim River	Republic of Serbia	43.393293 N	19.642978 E	9 Aug. 2016	Đuknić J.
MH549569	Simulium balcanicum/turgaicum 1	Pupa	Lim River	Republic of Serbia	43.393293 N	19.642978 E	9 Aug. 2016	Đuknić J.
MH549570	Simulium balcanicum/turgaicum 1	Pupa	Despotovica River	Republic of Serbia	43.975009 N	20.419295 E	10 Aug. 2016	Đuknić J.
MH549568	Simulium balcanicum/turgaicum 1	Larvae	Sava River	Republic of Croatia	45.759639 N	16.047861 E	4 Sept. 2015	Đuknić J.
MH587354	Simulium balcanicum/turgaicum 1	Larvae	Lim River	Republic of Serbia	43.393293 N	19.642978 E	9 Aug. 2016	Đuknić J.
MH549567	Simulium balcanicum/turgaicum 1	Pupa	Čemernica River	Republic of Serbia	43.909483 N	20.403543 E	10 Aug. 2016	Đuknić J.
MH587353	Simulium balcanicum/turgaicum 1	Pupa	Čemernica River	Republic of Serbia	43.909483 N	20.403543 E	10 Aug. 2016	Đuknić J.
MH587357	Simulium balcanicum	Larvae	Čemernica River	Republic of Serbia	43.909483 N	20.403543 E	10 Aug. 2016	Đuknić J.
MH587355	Simulium balcanicum	Pupa	Pek River	Republic of Serbia	44.490178 N	21.632081 E	20 Aug. 2015	Đuknić J.
MH513637	Simulium balcanicum	Pupa	Sava River	Republic of Slovenia	45.884078 N	15.640831 E	3 Sept. 2015	Đuknić J.
MH587358	Simulium balcanicum	Pupa	Semani River	Republic of Albania	40.750170 N	19.579350 E	26 Nov. 2016	Csányi B.
MH549561	Simulium balcanicum	Pupa	Sava River	Republic of Croatia	45.759639 N	16.047861 E	4 Sept. 2015	Đuknić J.
MH587359	Simulium balcanicum	Larvae	Despotovica River	Republic of Serbia	43.975009 N	20.419295 E	10 Aug. 2016	Đuknić J.
MH587360	Simulium balcanicum	Pupa	Marica River	Republic of Bulgaria	41.926976 N	25.933149 E	23 Sept. 2017	Đuknić J.
MH638296	Simulium balcanicum	Larvae	Semani River	Republic of Albania	40.750170 N	19.579350 E	26 Nov. 2016	Csányi B.
MH587356	Simulium balcanicum	Larvae	Sava River	Republic of Slovenia	45.884078 N	15.640831 E	3 Sept. 2015	Đuknić J.
MH587361	Simulium balcanicum	Pupa	Povelić River	Bosnia and Herzegovina	45.076826 N	17.491039 E	5 Oct. 2017	Đuknić J.
MH638295	Simulium balcanicum	Pupa	Isomotinis River	Hellenic Republic	41.052367 N	25.618511 E	23 Sept. 2017	Đuknić J.
MH549566	Simulium turgaicum 2	Larvae	Una River	Bosnia and Herzegovina	45.051033 N	16.379216 E	4 Oct. 2017	Đuknić J.
MH549565	Simulium turgaicum 2	Pupa	Sava River	Republic of Slovenia	45.884078 N	15.640831 E	3 Sept. 2015	Đuknić J.
MH549562	Simulium lineatum 2	Larvae	Sava River	Republic of Slovenia	45.884078 N	15.640831 E	3 Sept. 2015	Đuknić J.
MH549564	Simulium lineatum 2	Larvae	Sava River	Republic of Croatia	45.759639 N	16.047861 E	4 Sept. 2015	Đuknić J.
MH549563	Simulium lineatum 2	Pupa	Sava River	Republic of Croatia	45.759639 N	16.047861 E	4 Sept. 2015	Đuknić J.
MH151331	Simulium paraequinum 1	Pupa	Morača River	Montenegro	42.476927 N	19.304713 E	8 Aug. 2016	Đuknić J.
MH002239	Simulium paraequinum 1	Pupa	Čemernica River	Republic of Serbia	43.909483 N	20.403543 E	10 Aug. 2016	Đuknić J.
MH151332	Simulium paraequinum 1	Pupa	Lim River	Republic of Serbia	43.393293 N	19.642978 E	9 Aug. 2016	Đuknić J.
MH215257	Simulium paraequinum 1	Larvae	Morača River	Montenegro	42.476927 N	19.304713 E	8 Aug. 2016	Đuknić J.
MH549552	Simulium equinum	Pupa	Despotovica River	Republic of Serbia	43.975009 N	20.419295 E	10 Aug. 2016	Đuknić J.
MH215258	Simulium equinum	Pupa	Lim River	Republic of Serbia	43.393293 N	19.642978 E	9 Aug. 2016	Đuknić J.
MH549556	Simulium equinum	Larvae	Čemernica River	Republic of Serbia	43.909483 N	20.403543 E	10 Aug. 2016	Đuknić J.
MH549549	Simulium equinum	Pupa	Leva reka River	FYROM	41.142257 N	21.000276 E	22 Jun. 2017	Đuknić J.
MH549553	Simulium equinum	Pupa	Dragovištica River	Republic of Serbia	42.427797 N	22.520684 E	7 Oct. 2014	Đuknić J.
MH549551	Simulium equinum	Pupa	Čemernica River	Republic of Serbia	43.909483 N	20.403543 E	10 Aug. 2016	Đuknić J.
MH549554	Simulium equinum	Pupa	Dragovištica River	Republic of Serbia	42.427797 N	22.520684 E	7 Oct. 2014	Đuknić J.
MH549557	Simulium equinum	Pupa	Vrbas River	Bosnia and Herzegovina	44.624718 N	17.152975 E	5 Oct. 2017	Đuknić J.
MH549548	Simulium equinum	Larvae	Čemernica River	Republic of Serbia	43.909483 N	20.403543 E	10 Aug. 2016	Đuknić J.
MH549550	Simulium equinum	Pupa	Struma River	Republic of Bulgaria	41.943417 N	23.097965 E	17 Sept. 2017	Đuknić J.
MH549555	Simulium equinum	Pupa	Šemnica River	FYROM	41.065155 N	21.258429 E	23 Jun. 2017	Đuknić J.
MH587348	Simulium pseudequinum B	Pupa	Isomotinis River	Hellenic Republic	41.052367 N	25.618511 E	23 Sept. 2017	Đuknić J.
MH549558	Simulium pseudequinum B	Pupa	Tara River	Montenegro	42.863386 N	19.527027 E	8 Aug. 2016	Đuknić J.
MH549559	Simulium pseudequinum B	Pupa	Vrbica River	Republic of Bulgaria	41.569597 N	25.391039 E	23 Sept. 2017	Đuknić J.
MH587352	Simulium pseudequinum B	Pupa	Orahovica River	Republic of Croatia	45.550243 N	17.899489 E	25 Sept. 2016	Đuknić J.
MH587351	Simulium pseudequinum B	Pupa	Despotovica River	Republic of Serbia	43.975009 N	20.419295 E	10 Aug. 2016	Đuknić J.
MH549560	Simulium pseudequinum B	Larvae	Despotovica River	Republic of Serbia	43.975009 N	20.419295 E	10 Aug. 2016	Đuknić J.
MH587347	Simulium pseudequinum B	Larvae	Orahovica River	Republic of Croatia	45.550243 N	17.899489 E	25 Sept. 2016	Đuknić J.
MH587349	Simulium pseudequinum B	Larvae	Chlomo River	Hellenic Republic	40.587002 N	23.485925 E	18 Sept. 2017	Đuknić J.
MH587350	Simulium pseudequinum B	Larvae	Vrbica River	Republic of Bulgaria	41.569597 N	25.391039 E	23 Sept. 2017	Đuknić J.

Species names were given according to this study results.

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.

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 dH2O, 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.

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.

Results

The sequence length obtained from the studied 226 specimens ranged from 638 to 713 bp. The Tamura 3-parameter model with the gamma distribution of variation between the nucleotide positions (Tamura 1992) was shown to describe the substitution pattern the best (Table 2), as it had the lowest BIC score.

Table 2.

Five nucleotide substitution models that best fit the input data

Model	BIC	lnL
T92 + G + I	16449.48257	−5530.584276
T92 + G	16464.97209	−5544.276406
GTR + G + I	16507.1061	−5523.711842
GTR + G	16518.38881	−5535.300565
TN93 + G + I	16622.25964	−5599.130714

BIC, Bayesian information criterion (Schwarz 1978); GTR, General Time Reversible (Tavaré 1986); lnL, log-likelihood value. I = model with invariant sites; Γ = model with gamma distributed evolutionary rates among sites; T92 = Tamura 3-parameter (Tamura 1992); TN93 = Tamura-Nei (Tamura and Nei 1993) .

The topology of the phylogenetic tree of Wilhelmia was comprised of 15 monophyletic clades within the subgenus, which did not fully corroborate the morphological species identification (Fig. 2). In the GMYC analysis, the ML number of delimited species was also 15. These species are identical to the clades on the phylogenetic tree (Fig. 2). Following phylogenetic and GMYC analyses, we chose to respect the monophyly of all recovered clades (potential species), and based on the references, provided them with operational names.

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.

Morphologically defined species were present as monophyletic only in the case of Simulium sergenti Edwards, 1923, S. lineatum, and S. paraequinum, with the latter two comprised of two clades. All other represented Wilhelmia species were not seen as monophyletic. For instance, paraphyletic S. pseudequinum included four clades and polyphyletic S. equinum included two clades. Furthermore, some sequences were shown to arise from originally misidentified individuals. The names of the clades in the tree were either given in agreement with previous authors (e.g., Inci et al. 2017), or tentatively. The clades of S. pseudequinum from western Europe (the United Kingdom, France, and Spain) were named S. pseudequinum W1 and S. pseudequinum W2; the sequences from the Balkan Peninsula were named S. pseudequinum B, and the sequences from eastern Europe and Asia (Armenia, Turkey, and Pakistan) were called S. pseudequinum E.

Two main highly supported monophyletic branches (with BI >0.90) can be seen within the Wilhelmia subgenus (Fig. 2). One branch consisted of the clades S. equinum, S. pseudequinum B, S. pseudequinum E, S. pseudequinum W1, S. pseudequinum W2, and S. equinum 2. The second branch consisted of sergenti, paraequinum, and lineatum lineages. However, the positions of the sergenti and paraequinum lineages varied among the phylogeny reconstruction methods. The lineatum lineage was very diverse and included six clades: S. lineatum 1, S. lineatum 2, S. balcanicum, S. balcanicum/turgaicum 1, S. turgaicum 1, and S. turgaicum 2. The clade S. balcanicum/turgaicum 1 contained sequences from individuals identified morphologically as S. balcanicum (all pupae had bifurcate gills, as can be seen in Supp Fig. S3 [online only]), and sequences from individuals (morphologically and karyologically) identified as S. turgaicum (part of clade S. turgaicum 1, Inci et al. 2017).

Samples from the Balkan Peninsula occurred within seven clades. The previous checklist of Wilhelmia in the Balkans included S. angustifurca, S. balcanicum, S. equinum, S. lineatum, S. paraequinum, and S. pseudequinum (Adler and Crosskey, 2018). Samples of S. angustifurca were not available, so they were not included in this study. Significantly, our molecular results show the endemic nature of S. pseudequinum (clade S. pseudequinum B) and the presence of S. turgaicum (clade S. turgaicum 2) in the Balkans.

Nucleotide diversity within the monophyletic clades ranged from 2.5% within S. paraequinum 2, to 10.72% within S. equinum 2 (Table 3). The COI gene revealed that a highest haplotype diversity (1.000) was within species S. paraequinum 1, S. pseudequinum W2, and S. equinum 2, while the lowest (0.629) was within S. turgaicum 1. The highest number of haplotypes (32) was recorded for S. equinum (Table 3). Mismatch distribution testing of population expansion of the clades is shown in Supp Fig. S2 (online only). The frequency graphs of the pairwise differences between alleles indicated multimodal mismatch distribution. The negative values of Tajima’s D and Fu’s Fs (observed in most clades) indicate low nucleotide diversity but high haplotype diversity.

Table 3.

Nucleotide diversity calculations and tests of neutrality

Clades	n	h	S	Hd	Pi	Tajima’s D	Fu’s Fs
Simulium turgaicum 1	15	5	4/555	0.629 ± 0.125	0.00264 ± 0.00052	0.61628	−0.571
Simulium balcanicum/turgaicum 1	15	12	20/547	0.962 ± 0.040	0.00963 ± 0.00105	−0.75871	−4.539
Simulium balcanicum	29	15	24/538	0.872 ± 0.053	0.00783 ± 0.00119	−1.21092	−4.569
Simulium turgaicum 2	9	6	13/612	0.889 ± 0.091	0.00654 ± 0.00136	−0.78119	−0.5
Simulium lineatum 2	13	8	10/531	0.859 ± 0.089	0.00377 ± 0.00075	−1.51552	−3.645
Simulium lineatum 1	11	8	11/552	0.891 ± 0.092	0.00580 ± 0.00096	−0.63656	−2.782
Simulium paraequinum 1	6	6	13/646	1.000 ± 0.096	0.00846 ± 0.00161	−0.24351	−1.987
Simulium paraequinum 2	10	5	5/525	0.756 ± 0.130	0.00250 ± 0.00068	−1.03527	−1.587
Simulium sergenti	15	7	10/658	0.781 ± 0.102	0.00324 ± 0.00088	−1.16596	−1.686
Simulium equinum	58	32	34/551	0.887 ± 0.038	0.00717 ± 0.00061	−1.62754	−24.298*
Simulium pseudequinum B	9	4	9/565	0.750 ± 0.112	0.00531 ± 0.00127	−0.43284	1.293
Simulium pseudequinum E	8	7	17/615	0.964 ± 0.077	0.00952 ± 0.00149	−0.55019	−1.547
Simulium pseudequinum W1	6	4	8/658	0.800 ± 0.172	0.00638 ± 0.00128	1.1717	0.851
Simulium pseudequinum W2	3	3	4/634	1.000 ± 0.272	0.00421 ± 0.00140	–	−0.341
Simulium equinum 2	9	9	23/632	1.000 ± 0.052	0.01072 ± 0.00146	−1.05159	−3.932

h = number of haplotypes; Hd = haplotype diversity ± SD; n = number of sequences; Pi = nucleotide diversity ± SD; S = number of segregating sites; Tajima’s D test and Fu’s Fs test.

*P < 0.02.

The inter-clade divergence for the COI sequence of Wilhelmia ranged from 1.78% (S. lineatum 1 vs. S. lineatum 2) to 18.93% (S. equinum 1 vs. S. paraequinum 1) (Table 4). All clades from the lineatum lineage (S. lineatum 1, S. lineatum 2, S. balcanicum, S. balcanicum/turgaicum 1, S. turgaicum 1, and S. turgaicum 2) displayed low genetic distances from each other (1.78–3.30%) (Table 4).

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)

Clades	1.	2.	3.	4.	5.	6.	7.	8.	9.	10.	11.	12.	13.	14.
1. Simulium turgaicum 1
2. Simulium balcanicum/ turgaicum1	0.0212
3. Simulium balcanicum	0.0250	0.0204
4. Simulium turgaicum 2	0.0308	0.0256	0.0288
5. Simulium lineatum 2	0.0330	0.0271	0.0303	0.0254
6. Simulium lineatum 1	0.0302	0.0275	0.0279	0.0262	0.0178
7. Simulium paraequinum 1	0.1307	0.1338	0.1329	0.1313	0.1367	0.1284
8. Simulium paraequinum 2	0.1142	0.1202	0.1216	0.1181	0.1159	0.1152	0.0510
9. Simulium sergenti	0.1273	0.1268	0.1243	0.1241	0.1224	0.1212	0.1487	0.1293
10. Simulium equinum 1	0.1624	0.1711	0.1668	0.1656	0.1509	0.1551	0.1893	0.1756	0.1613
11. Simulium pseudequinum B	0.1554	0.1565	0.1512	0.1455	0.1492	0.1667	0.2031	0.1734	0.1552	0.0862
12. Simulium pseudequinum E	0.1487	0.1486	0.1425	0.1459	0.1333	0.1350	0.1550	0.1550	0.1537	0.0842	0.0833
13. Simulium pseudequinum W1	0.1490	0.1545	0.1505	0.1540	0.1481	0.1475	0.1783	0.1690	0.1345	0.0993	0.0993	0.0990
14. Simulium pseudequinum W2	0.1355	0.1486	0.1518	0.1463	0.1413	0.1373	0.1741	0.1552	0.1360	0.0937	0.1071	0.0898	0.0259
15. Simulium equinum 2	0.1364	0.1507	0.1482	0.1358	0.1406	0.1346	0.1618	0.1555	0.1610	0.1779	0.1572	0.1399	0.1404	0.1203

A total of 122 haplotypes of Wilhelmia species was recognized in DnaSP (Fig. 3). After applying the star contraction method, the number of haplotypes was reduced to 80. The minimal number of mutations (5–7) was recorded among the clades of lineatum lineage, while the maximal number of mutations was found between the equinum and lineatum branches (54), S. equinum 2 and the other clades from the equinum branch (42), and between sergenti and lineatum lineages (41).

Fig. 3.

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

According to the time-event analyses, the evolution of the subgenus Wilhelmia started in the Late Cretaceous (111–67 Ma) (Fig. 4). The first diversification within the subgenus into two branches started by the end of the Cretaceous and during the Paleocene (76–46 Ma) (Figs. 4 and 5A). Many distant islands in Tethys (distances up to 1,500 km) surrounded by continents, as well as a variety of climatic zones characterized the western Palearctic at that time and provided a great potential for allopatric speciation events. From the potential islands of origin, branches and lineages of species could spread their distribution in a stepwise manner. Possible allopatric speciation can also be observed within branches. According to our results, in morphologically uniform S. pseudequinum, diversification among geographically distant clades occurred in 46–21 Ma during the Oligocene (Fig. 5B).

Fig. 4.

Bayesian phylogenetic tree based on the COI gene of Wilhelmia clades, with S. notiale, Culicoides anophelis, and Thaumalea testacea as the outgroup. Numbers above the branches represent the mean time value for the branching event, while the numbers below the branches represent a confidence range of the branching time.

Fig. 5.

Paleogeological maps showing the location and environment/climate of land masses (following Kazmin and Natapov 1998). A) Maastrichtian age (~70 Ma), time of first branching within the subgenus; B) Oligocene (~30Ma), time of diversification among S. pseudequinum.

Discussion

This study applied molecular barcoding in addition to morphological identification of immature life stages of black flies in the Balkan Peninsula. Samples of individuals collected in the Balkans were analyzed together with all available Wilhelmia sequences within the phylogenetic frame of the Simulium subgenus Wilhelmia. In this study, 15 monophyletic clades were recognized (potential species) out of 7 morphologically identified species.

Our phylogenetic analyses discovered two major branches within the subgenus Wilhelmia. One branch included the morphospecies S. equinum and S. pseudequinum, while the other included morphospecies S. sergenti, S. paraequinum, S. lineatum, S. turgaicum, and S. balcanicum. Identical divergence among these species assemblages was previously recorded for larval polytene chromosome data (Weber and Grunewald 1989, Petrova et al. 2003, Chubareva et al. 2007, Huang et al. 2012, Adler et al. 2015), which lead to the naming of the assemblages (Wilhelmia) equina group and (Wilhelmia) salopiensis group, respectively (Chubareva et al. 2007). As Wilhelmia salopiense Edwards, 1927 is regarded as a synonym of S. lineatum (Adler and Crosskey 2018), we prefer to call this branch the lineatum branch. The positions of some lineages within the lineatum branch differed among ML, MP, and Bayesian phylogenetic trees. We explained the evolution of Wilhelmia following the Bayesian tree topology, as it was in accordance with previous evolutionary relations obtained by comparing chromosomal characters.

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.

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.

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).

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.

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.

Biodiversity in Balkan Peninsula

The inventory given by Adler and Crosskey (2018) contains six species of Simulium subg. Wilhelmia for the Balkans (60% of species present in Europe): S. balcanicum, S. lineatum, S. paraequinum, S. equinum, S. pseudequinum, and S. angustifurca (the last one was not present in this study). The Balkan Peninsula is usually seen as one of southern Europe’s biodiversity hotspots. However, based on the current inventory, the Balkan Peninsula does not bear extraordinary Wilhelmia species richness in comparison with the surrounding regions (central Europe—six, Anatolia and Caucasus—seven; Adler and Crosskey 2018). Since significant cryptic diversity within the Wilhelmia subgenus is continuously being discovered, the taxa richness in the mentioned area could be quite different.

Our sampling and barcoding of Wilhelmia revealed the presence of seven clades (potential species) in the studied Balkan area. The morphospecies S. equinum, S. pseudequinum, S. paraequinum, S. lineatum, and S. turgaicum were each represented by one of the clades (S. equinum, S. pseudequinum B, S. paraequinum 1, S. lineatum 2, and S. turgaicum 2). However, samples morphologically identified as S. balcanicum belonged to two clades: S. balcanicum and S. balcanicum/turgaicum 1. Clades S. balcanicum, S. pseudequinum B, and S. equinum were the most widely present in the Balkans.

Three morphospecies of the lineatum lineage (S. balcanicum, S. lineatum, and S. turgaicum) differ in their distributions, according to Adler et al. (2015). Widespread species S. balcanicum increase in prevalence eastward. The species S. lineatum is present throughout Europe and its eastern areal boundary is in the Balkans, in Bulgaria, while it has been taken over by S. turgaicum in the Middle East. Adler et al. (2015) assumed that there is possible overlap of S. lineatum and S. turgaicum distribution, potentially somewhere in southern Russia or in the Balkans, perhaps Bulgaria. The same authors have remarked on some of the ecological differences between these taxa (S. lineatum is a lowland species, present below 500 m and S. turgaicum occurs above 900 m). Also, when S. balcanicum and S. lineatum overlap geographically, they were reported to typically occur in different rivers (Adler et al. 2015). In our study, S. turgaicum (sequences from the clade S. turgaicum 2) was recorded in the Balkans for the first time. New sequences of Wilhelmia species from the Balkan Peninsula showed that all three species (S. balcanicum, S. lineatum, and S. turgaicum) can be found in sympatry, which was confirmed by samples from the Sava River in Slovenia (elevation 132 m). These findings showed that the western boundary of S. turgaicum distribution extended far more to the west, while the western part of the Balkan Peninsula was indeed seen as the southeastern boundary of the range of S. lineatum.

The most intriguing result of this study was the geographical differentiation in paraphyletic morphospecies S. pseudequinum. Because the divergence of its clades took place during the Oligocene, each of them could be endemic to its area of origin. Thus, if it proves to be correct, S. pseudequinum B would be endemic to the Balkan Peninsula. Further molecular studies at areas of potential contact between S. pseudequinum clades are needed to better interpret their distributions and evolutionary histories.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.