## Conclusions

This study provides clarification of the distribution pattern of species belonging to the Obsoletus/Scoticus Complex, using a dataset based on samples from western Palaearctic and Mediterranean transect. Strong variations in latitudinal cryptic species diversity was observed. This study clarifies the phylogenetic relationships between species belonging to the Obsoletus/Scoticus Complex. We identified and validated five MOTUs, *C. obsoletus*, *C. scoticus*, *C. obsoletus* clade O2, *C. obsoletus* clade dark and a MOTUs corresponding to an unidentified species. The latter three species have not been formally described but our results confirm that they should be considered as species in their own right. More detailed studies of their morphology and ecology are needed to provide more detailed descriptions of these species. Furthermore, our results raise questions concerning the taxonomic status of *C. montanus*, which was previously considered as a taxonomically valid species.

## Supplementary information

**Supplementary information** accompanies this paper at https://doi.org/10.1186/s13071-020-04114-1.

**Additional file 1: Table S1.** Information on adult female *Culicoides* sampling sites and results of specific assignation using *cox1*.

**Additional file 2: Table S2** Primers used for PCRs and sequencing in this study.

**Additional file 3: Table S3.** Reference sequences used for specific assignation.

**Additional file 4: Table S4.** Information on all *Culicoides* sequenced and GenBank accession numbers.

**Additional file 5: Figure S1.** Maximum likelihood phylogenetic tree using 16S rDNA. Values at the nodes represent bootstrap values greater than 800 for the most important nodes (1000 replicates).

**Additional file 6: Figure S2.** Maximum likelihood phylogenetic tree using *cox1*. Values at the nodes represent bootstrap values greater than 800 for the most important nodes (1000 replicates).

**Additional file 7: Figure S3.** Maximum likelihood phylogenetic tree using 28S rDNA. Values at the nodes represent bootstrap values greater than 800 for the most important nodes (1000 replicates).

## Abbreviations

BTV: bluetongue virus; SBV: Schmallenberg virus; *cytb*: cytochrome oxidase b; ITS1: internal transcribed spacer region 1; ITS2: internal transcribed spacer region 2; *cox1*: cytochrome c oxidase subunit 1; N: number of individuals; h: number of haplotypes; Hd: haplotype (gene) diversity; nt: nucleotide diversity; nd: not determined; S: number of polymorphic sites; bGMYC: Bayesian General Mixed Yule Coalescent; MOTUs: molecular operational taxonomic units; bPTP: Bayesian Poisson Tree Processes.

## Acknowledgements

The authors are grateful to all European partners who assisted in *Culicoides* sampling and shipping. We acknowledge Bruno Mathieu for his friendly cooperation and for allowing us access to use some of his sequences as references for our study.

## Authors' contributions

AM, KH and CG designed the study. AJ, AD, AC, BVP, DWR, DS, DW, DPu, DPe, EV, FJ, HK, IPF, JL, JN, JMP, JS, KRS, KK, CLC, ML, MBo, MG, MBi, ME, MR, MO,

MAMC, RB, REP, SC, ST, SB, SS, SMS, ZO, ZS and ZV provided samples of *Culicoides* collected in their respective countries. IR, MD contributed the identification of *Culicoides*. AM, LG, LP and LT performed DNA extraction, amplification, quality control and sequencing. BM provided some sequences as references for our study. AM analysed the data. KH and CG contributed to the manuscript written first by AM. All authors read and approved on the final manuscript.

## Funding

This study was partially funded by the EU grant H2020-727393 PALE-Blu and by the VectorNet project (OC/EFSA/AHAW/2013/02-FWC1) funded by the European Centre for Disease Prevention and Control (ECDC) and the European Food Safety Authority (EFSA). The contents of this publication are the sole responsibility of the authors and do not necessarily reflect the views of the European Commission, of the ECDC or of the EFSA.

## Availability of data and materials

All data generated or analysed during this study are included in this published article and its additional files. The newly generated sequences were submitted in the GenBank database under the accession numbers MT170026-MT173788.

## Ethics approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Competing interests

The authors declare that they have no competing interests.

## Author details

¹ ASTRE, Univ Montpellier, Cirad, INRAE, Montpellier, France. ² Cirad, UMR ASTRE, 34398 Montpellier, France. ³ Cirad, UMR ASTRE, 97170 Petit-Bourg, Guadeloupe, France. ⁴ Institut Agronomique et Vétérinaire Hassan II, Unité Parasitologie et Maladies Parasitaires, 10100 Rabat, Morocco. ⁵ Institute of Parasitology and Tropical Pathology of Strasbourg, Université de Strasbourg, DIHP UR 7292, 67000 Strasbourg, France. ⁶ Centre for Monitoring of Vectors, National Reference Centre, Netherlands Food and Consumer Product Safety Authority, Wageningen, The Netherlands. ⁷ Veterinary Control Central Research Institute, Ankara, Turkey. ⁸ Department of Parasitology and Parasitic Diseases, Faculty of Veterinary Medicine, Ss. Cyril and Methodius University in Skopje, Skopje, Republic of North Macedonia. ⁹ Centre for Ecology, Centre for Ecology & Hydrology, Wallingford OX10 8BB, UK. ¹⁰ CIISA-Centro de Investigação Interdisciplinar em Sanidade Animal, Faculdade de Medicina Veterinária, Universidade de Lisboa, Avenida da Universidade Técnica, 1300-477 Lisboa, Portugal. ¹¹ Department of Parasitology-Parasitic Diseases, Entomology & Bee Health, Veterinary Centre of Athens, Athens, Greece.

¹² Leibniz-Centre for Agricultural Landscape Research, Müncheberg, Germany.

¹³ Faculty of Agriculture, University of Novi Sad, Novi Sad, Serbia. ¹⁴ National Centre for Vector Entomology, Institute of Parasitology, University of Zürich, Zürich, Switzerland. ¹⁵ Friedrich-Loeffler-Institut, Federal Research Institute for Animal Health, Greifswald, Germany. ¹⁶ Department of Animal Pathology, AgriFood Institute of Aragón (IA2) Veterinary Faculty, 50013 Zaragoza, Spain.

¹⁷ Departamento de Microbiología, Laboratorio de Producción y Sanidad Animal de Granada, Junta de Andalucía, Granada, Spain. ¹⁸ Doñana Biological Station, CSIC, Sevilla, Spain. ¹⁹ Centro de Investigación Biomédica en Red de Epidemiología y Salud Pública (CIBERESP), Madrid, Spain. ²⁰ Centre for Ecology & Hydrology, Edinburgh OX10 8BB, UK. ²¹ Department of Virology, University of Helsinki, Medicum, Haartmaninkatu 3, Helsinki 00014, Finland. ²² National Veterinary Research Institute, Pulawy, Poland. ²³ Istituto Zooprofilattico Sperimentale dell'Abruzzo e del Molise 'G. Caporale', Campo Boario, 64100 Teramo, Italy. ²⁴ The Pirbright Institute, Pirbright, UK. ²⁵ Department of Epidemiology and Population Health, Institute of Infection and Global Health, University of Liverpool, Leahurst, Chester High Road, Neston, Cheshire CH64 7TE, UK.

²⁶ Applied Zoology and Animal Conservation Research Group, University of the Balearic Islands UIB, Palma, Spain. ²⁷ University of Copenhagen, Copenhagen, Denmark. ²⁸ National Diagnostic and Research Veterinary Medical Institute, Sofia, Bulgaria. ²⁹ Norwegian Veterinary Institute, Oslo, Norway.

³⁰ Institute of Food safety, Animal Health and Environment 'BIOR', Riga, Latvia.
Received: 13 January 2020 Accepted: 29 April 2020

Published online: 20 May 2020

References

1. Carpenter S, Groschup MH, Garros C, Felippe-Bauer ML, Purse BV. *Culicoides* biting midges, arboviruses and public health in Europe. Antiviral Res. 2013;100:102-13.

2. Wilson AJ, Mellor PS. Bluetongue in Europe: past, present and future. Phil Trans R Soc Lond B. 2009;364:2669-81.

3. Carpenter S, Wilson A, Mellor PS. *Culicoides* and the emergence of bluetongue virus in northern Europe. Trends Microbiol. 2009;17:172-8.

4. Hoffmann B, Eschbaumer M, Beer M. Real-time quantitative reverse transcription-PCR assays specifically detecting bluetongue virus serotypes 1, 6, and 8. J Clin Microbiol. 2009;47:2992-4.

5. Savini G, MacLachlan NJ, Sanchez-Vizcaino JM, Zientara S. Vaccines against bluetongue in Europe. Comp Immunol Microbiol Infect Dis. 2008;31:101-20.

6. Carpenter S, Lunt HL, Arav D, Venter GJ, Mellor PS. Oral susceptibility to bluetongue virus of *Culicoides* (Diptera: Ceratopogonidae) from the United Kingdom. J Med Entomol. 2006;43:73-8.

7. Carpenter S, Mcarthur C, Selby R, Ward R, Nolan DV, Luntz AJ, et al. Experimental infection studies of UK *Culicoides* species midges with bluetongue virus serotypes 8 and 9. Vet Rec. 2008;163:589-92.

8. Mathieu B. Les espèces de *Culicoides* du sous-genre *Avaritia* (Diptera: Ceratopogonidae) dans le monde: révision systématique et taxonomique des espèces d'intérêt dans la transmission d'Orbivirus. École doctorale Sciences de la vie et de la santé Université de Strasbourg; 2011.

9. Meiswinkel R, Gomulski LM, Delecolle JC, Goffredo M, Gasperi G. The taxonomy of *Culicoides* vector complexes—unfinished business. Vet Ital. 2004;40:151-9.

10. Hajd Henni L, Sauvage F, Ninio C, Depaquit J, Augot D. Wing geometry as a tool for discrimination of Obsoletus group (Diptera: Ceratopogonidae: *Culicoides*) in France. Infect Genet Evol. 2014;21:110-7.

11. Mehlhorn H, Walldorf V, Klimpel S, Jahn B, Jaeger F, Eschweiler J, et al. First occurrence of *Culicoides obsoletus*-transmitted bluetongue virus epidemic in central Europe. Parasitol Res. 2007;101:219-28.

12. Savini G, Goffredo M, Monaco F, Di Gennaro A, Cafiero MA, Baldi L, et al. Bluetongue virus isolations from midges belonging to the *Obsoletus* complex (*Culicoides*, Diptera: Ceratopogonidae) in Italy. Vet Rec. 2005;157:133-9.

13. Borkent A. World species of biting midges (Diptera: Ceratopogonidae). Salmon Arm: American Museum of Natural History, and Instituto Nacional de Biodiversidad; 2016.

14. Gomulski LM, Meiswinkel R, Delecolle JC, Goffredo M, Gasperi G. Phylogenetic relationships of the subgenus *Avaritia* Fox, 1955 including *Culicoides obsoletus* (Diptera, Ceratopogonidae) in Italy based on internal transcribed spacer 2 ribosomal DNA sequences. Syst Entomol. 2005;30:619-31.

15. Ander M, Troell K, Chirico J. Barcoding of biting midges in the genus *Culicoides*: a tool for species determination. Med Vet Entomol. 2013;27:323-31.

16. Kluiters G, Pages N, Carpenter S, Gardes L, Guis H, Baylis M, et al. Morphometric discrimination of two sympatric sibling species in the Palaearctic region, *Culicoides obsoletus* Meigen and *C. scoticus* Downes & Kettle (Diptera: Ceratopogonidae), vectors of bluetongue and Schmallenberg viruses. Parasit Vectors. 2016;9:262.

17. Ninio C, Augot D, Delecolle JC, Dufour B, Depaquit J. Contribution to the knowledge of *Culicoides* (Diptera: Ceratopogonidae) host preferences in France. Parasitol Res. 2011;108:657-63.

18. Martinez-de la Puente J, Figuerola J, Soriguer R. Fur or feather? Feeding preferences of species of *Culicoides* biting midges in Europe. Trends Parasitol. 2015;31:16-22.

19. Garros C, Gardes L, Allene X, Rakotoarivony I, Viennet E, Rossi S, et al. Adaptation of a species-specific multiplex PCR assay for the identification of blood meal source in *Culicoides* (Ceratopogonidae: Diptera): applications on Palaearctic biting midge species, vectors of Orbiviruses. Infect Genet Evol. 2011;11:1103-10.

20. Kettle DS, Lawson JWH. The early stages of British biting midges *Culicoides* Latreille (Diptera: Ceratopogonidae) and allied Genera. Bull Entomol Res. 2009;43:421-67.

21. Meiswinkel R, De Bree F, Bossers-De Vries R, Elbers AR. An unrecognized species of the *Culicoides obsoletus* complex feeding on livestock in The Netherlands. Vet Parasitol. 2015;207:324-8.

22. Elbers AR, Meiswinkel R. *Culicoides* (Diptera: Ceratopogonidae) host preferences and biting rates in the Netherlands: comparing cattle, sheep and the black-light suction trap. Vet Parasitol. 2014;205:330-7.

23. Harbach RE. Review of the internal classification of the genus *Anopheles* (Diptera, Culicidae) - the foundation for comparative systematics and phylogenetic research. Bull Entomol Res. 1994;84:331-42.

24. Debeila TJ. Characterisatioin of selected *Culicoides* population in South Africa using genetic markers. PhD Thesis, Faculty of Veterinary Science, University of Pretoria; 2010.

25. Tabachnick W. *Culicoides* and the global epidemiology of bluetongue virus infection. Vet Ital. 2004;40:144-50.

26. Bakhoum MT, Labuschagne K, Huber K, Fall M, Mathieu B, Venter G, et al. Phylogenetic relationships and molecular delimitation of *Culicoides* Latreille (Diptera: Ceratopogonidae) species in the Afrotropical region: interest for the subgenus *Avaritia*. Syst Entomol. 2018;43:355-71.

27. Cetre-Sossah C, Baldet T, Delecolle JC, Mathieu B, Perrin A, Grillet C, et al. Molecular detection of *Culicoides* spp. and *Culicoides imicola*, the principal vector of bluetongue (BT) and African horse sickness (AHS) in Africa and Europe. Vet Res. 2004;35:325-37.

28. Gomulski LM, Meiswinkel R, Delecolle JC, Goffredo M, Gasperi G. Phylogeny of the subgenus *Culicoides* and related species in Italy, inferred from internal transcribed spacer 2 ribosomal DNA sequences. Med Vet Entomol. 2006;20:229-38.

29. Nolan DV, Carpenter S, Barber J, Mellor PS, Dallas JF, Mordue Luntz AJ, et al. Rapid diagnostic PCR assays for members of the *Culicoides obsoletus* and *Culicoides pulicaris* species complexes, implicated vectors of bluetongue virus in Europe. Vet Microbiol. 2007;124:82-94.

30. Augot D, Sauvage F, Jouet D, Simphal E, Veuille M, Couloux A, et al. Discrimination of *Culicoides obsoletus* and *Culicoides scoticus*, potential bluetongue vectors, by morphometrical and mitochondrial cytochrome oxidase subunit I analysis. Infect Genet Evol. 2010;10:629-37.

31. Pages N, Munoz-Munoz F, Talavera S, Sarto V, Lorca C, Nunez JI. Identification of cryptic species of *Culicoides* (Diptera: Ceratopogonidae) in the subgenus *Culicoides* and development of species-specific PCR assays based on barcode regions. Vet Parasitol. 2009;165:298-310.

32. Schwenkenbecher JM, Mordue AJ, Piertney SB. Phylogenetic analysis indicates that *Culicoides dewulff* should not be considered part of the *Culicoides obsoletus* complex. Bull Entomol Res. 2009;99:371-5.

33. Sarto i Monteys V, Ventura D, Pages N, Aranda C, Escosa R. Expansion of *Culicoides imicola*, the main bluetongue virus vector in Europe, into Catalonia, Spain. Vet Rec. 2005;156:415-7.

34. Harrup LE, Bellis GA, Balenghien T, Garros C. *Culicoides* Latreille (Diptera: Ceratopogonidae) taxonomy: current challenges and future directions. Infect Genet Evol. 2015;30:249-66.

35. Delecolle JC. Nouvelle contribution à l'étude systématique et iconographique des espèces du genre *Culicoides* (Diptéra: Ceratopogonidae) du Nord-Est de la France. PhD Thesis, University of Strasbourg, Strasbourg; 1985.

36. Campbell JA, Pelham-Ch Clinton EC. Taxonomic review of the British species of *Culicoides* Latreille (Diptera, Ceratopogonidae). Proc R Entomol Soc. 1960;67:181-302.

37. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792-7.

38. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451-2.

39. Castresana J. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol Biol Evol. 2000;17:540-52.

40. Wickham H. ggplot2: Elegant graphics for data analysis. J R Stat Soc A Stat. 2009;174:245-50.

41. Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol. 1980;16:111-20.

42. Nei M, Rooney AP. Concerted and birth-and-death evolution of multi-gene families. Annu Rev Genet. 2005;39:121-52.
43. Kumar S, Stecher G, Li M, Knyaz C, Tarnura K. MEGA X: molecular evolutionary genetics analysis across computing platforms. Mol Biol Evol. 2018;35:1547-9.

44. Reid NM, Carstens BC. Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yule-coalescent model. BMC Evol Biol. 2012;12:196.

45. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969-73.

46. Fujisawa T, Barraclough TG. Delimiting species using single-locus data and the Generalized Mixed Yule Coalescent approach: a revised method and evaluation on simulated data sets. Syst Biol. 2013;62:707-24.

47. Zhang J, Kapli P, Pavlidis P, Stamatakis A. A general species delimitation method with applications to phylogenetic placements. Bioinformatics. 2013;29:2869-76.

48. Crozier RH, Crozier YC. The mitochondrial genome of the honeybee *Apis mellifera*: complete sequence and genome organization. Genetics. 1993;133:97-117.

49. Möhlmann TWR, Bekendam AM, van Kemenade I, Wernergren U, Favia G, Takken W, et al. Latitudinal diversity of biting midge species within the Obsoletus group across three habitats in Europe. Med Vet Entomol. 2019;33:420-6.

50. Goffredo M, Meiswinkel R, Federici V, Di Nicola F, Mancini G, Ippoliti C, et al. The "*Culicoides obsoletus* group" in Italy: relative abundance, geographic range, and role as vector for bluetongue virus. Vet Ital. 2016;52:235-41.

51. Ippoliti C, Goffredo M, Salini R, Bocci M, Quaglia M, Federici V, et al. Ecolimatic indicators for three *Culicoides* species of the *Obsoletus* complex in Italy. Vet Ital. 2016;52:213-22.

52. Ander M, Meiswinkel R, Chirico J. Seasonal dynamics of biting midges (Diptera: Ceratopogonidae: *Culicoides*), the potential vectors of blue-tongue virus, in Sweden. Vet Parasitol. 2012;184:59-67.

53. Allen AP, Gillooly JF. Assessing latitudinal gradients in speciation rates and biodiversity at the global scale. Ecol Lett. 2006;9:947-54.

54. Wenk CE, Kaufmann C, Schaffner F, Mathis A. Molecular characterization of Swiss Ceratopogonidae (Diptera) and evaluation of real-time PCR assays for the identification of *Culicoides* biting midges. Vet Parasitol. 2012;184:258-66.

55. Augot D, Mathieu B, Hadj-Henni L, Barriel V, Zapata Mena S, Smolis S, et al. Molecular phylogeny of 42 species of *Culicoides* (Diptera, Ceratopogonidae) from three continents. Parasite. 2017;24:23.

56. Kirkeby C, Dominik P. *Culicoides (Avaritia) gornostaevae* Mirzaeva, 1984 (Diptera: Ceratopogonidae)—a possible vector species of the *Obsoletus* group new to the European fauna. Parasit Vectors. 2014;7:445.

57. Mathieu B, Delecolle JC, Garros C, Balenghien T, Setier-Rio ML, Candolfi E, et al. Simultaneous quantification of the relative abundance of species complex members: application to *Culicoides obsoletus* and *Culicoides scoticus* (Diptera: Ceratopogonidae), potential vectors of bluetongue virus. Vet Parasitol. 2011;182:297-306.

58. Pages N, Sarto IVM. Differentiation of *Culicoides obsoletus* and *Culicoides scoticus* (Diptera: Ceratopogonidae) based on mitochondrial cytochrome oxidase subunit I. J Med Entomol. 2005;42:1026-34.

59. Kiehl E, Walldorf V, Klimpel S, Al-Quraishy S, Mehlhorn H. The European vectors of Bluetongue virus: are there species complexes, single species or races in *Culicoides obsoletus* and *C. pulicaris* detectable by sequencing ITS-1, ITS-2 and 18S-rDNA? Parasitol Res. 2009;105:331-6.

60. Esselstyn JA, Evans BJ, Sedlock JL, Anwarali Khan FA, Heaney LR. Single-locus species delimitation: a test of the mixed Yule-coalescent model, with an empirical application to Philippine round-leaf bats. Proc R Entomol Soc. 2012;279:3678-86.

61. Talavera G, Dinca V, Vila R. Factors affecting species delimitations with the GMYC model: insights from a butterfly survey. Methods Ecol Evol 2013;4:1101-10.

62. Dellicour S, Flot JF. Delimiting species-poor data sets using single molecular markers: a study of barcode gaps, haplowebs and GMYC. Syst Biol. 2015;64:900-8.

63. Rosenberg NA. The probability of topological concordance of gene trees and species trees. Theor Popul Biol. 2002;61:225-47.

64. Nichols R. Gene trees and species trees are not the same. Trends Ecol Evol. 2001;16:358-64.

65. Pentinsaari M, Vos R, Mutanen M. Algorithmic single-locus species delimitation: effects of sampling effort, variation and nonmonophyly in four methods and 1870 species of beetles. Mol Ecol Resour. 2017;17:393-404.

66. Fujita MK, Leache AD, Burbrink FT, McGuire JA, Moritz C. Coalescent-based species delimitation in an integrative taxonomy. Trends Ecol Evol. 2012;27:480-8.

67. Leache AD, Fujita MK. Bayesian species delimitation in West African forest geckos (*Hemidactylus fasciatus*). Proc R Entomol Soc. 2010;277:3071-7.

68. Boykin LM, Schutze MK, Krosch MN, Chomič A, Chapman TA, Englezou A, et al. Multi-gene phylogenetic analysis of south-east Asian pest members of the *Bactrocera dorsalis* species complex (Diptera: Tephritidae) does not support current taxonomy. J Appl Entomol. 2014;138:235-53.

69. Dinca V, Lukhtanov VA, Talavera G, Vila R. Unexpected layers of cryptic diversity in wood white *Leptidea butterflies*. Nat Commun. 2011;2:324.

70. Hsieh CH, Ko CC, Chung CH, Wang HY. Multilocus approach to clarify species status and the divergence history of the *Bemisia tabaci* (Hemiptera: Aleyrodidae) species complex. Mol Phylogenet Evol. 2014;76:172-80.

71. Song J-H, Ahn K-J. Species delimitation in the *Aleochara fucicola* species complex (Coleoptera: Staphylinidae: Aleocharinae) and its phylogenetic relationships. Zool Scr. 2014;43:629-40.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_30.png" /></figure>
Supplementary information

<table><caption>**Tableau S1.** Information on adult female *Culicoides* sampling sites and results of specific assignation using cox1.</caption><thead><tr><th>Locations names</th><th>Country</th><th>Latitude</th><th>Longitude</th><th>Year of collection</th><th>Number of Culicoides obsoletus</th><th>Number of Culicoides montanus</th><th>Number of Culicoides obsoletus clade O2</th><th>Number of Culicoides scoticus</th><th>Number of Culicoides obsoletus clade dark</th><th>Number of Culicoides scoticus clade 2</th><th>Number of Culicoides sp.</th></tr></thead><tbody><tr><td>Montana</td><td>Bulgaria</td><td>43,41</td><td>23,23</td><td>2016</td><td>146</td><td>0</td><td>0</td><td>8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Nikolovo</td><td>Bulgaria</td><td>43,87</td><td>26,12</td><td>2016</td><td>146</td><td>2</td><td>0</td><td>15</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Dyssevej</td><td>Denmark</td><td>55,20</td><td>12,11</td><td>2015</td><td>110</td><td>0</td><td>21</td><td>6</td><td>1</td><td>0</td><td>0</td></tr><tr><td>Groddevej</td><td>Denmark</td><td>56,13</td><td>9,18</td><td>2015</td><td>43</td><td>0</td><td>0</td><td>10</td><td>0</td><td>1</td><td>0</td></tr><tr><td>Lægdsgårdsvej</td><td>Denmark</td><td>56,12</td><td>9,12</td><td>2015</td><td>89</td><td>0</td><td>30</td><td>20</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Leestrupvej</td><td>Denmark</td><td>55,22</td><td>12,05</td><td>2015</td><td>91</td><td>0</td><td>4</td><td>85</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Inari</td><td>Finland</td><td>68,67</td><td>27,56</td><td>2015</td><td>8</td><td>0</td><td>49</td><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>Pelkosennemi</td><td>Finland</td><td>67,18</td><td>27,46</td><td>2015</td><td>0</td><td>0</td><td>90</td><td>6</td><td>0</td><td>15</td><td>0</td></tr><tr><td>Sotkamo</td><td>Finland</td><td>64,04</td><td>28,72</td><td>2015</td><td>19</td><td>0</td><td>8</td><td>36</td><td>4</td><td>0</td><td>0</td></tr><tr><td>Barjols</td><td>France</td><td>43,54</td><td>5,98</td><td>2016</td><td>13</td><td>0</td><td>0</td><td>14</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Bonnay</td><td>France</td><td>49,93</td><td>2,51</td><td>2016</td><td>39</td><td>0</td><td>0</td><td>6</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Les Hitres</td><td>France</td><td>45,75</td><td>4,07</td><td>2016</td><td>43</td><td>0</td><td>0</td><td>4</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Reygnac</td><td>France</td><td>45,34</td><td>1,82</td><td>2016</td><td>11</td><td>0</td><td>7</td><td>0</td><td>12</td><td>0</td><td>0</td></tr><tr><td>Sucé-sur-Erdre</td><td>France</td><td>47,32</td><td>-1,54</td><td>2016</td><td>8</td><td>0</td><td>17</td><td>15</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Peetzig</td><td>Germany</td><td>53,08</td><td>13,94</td><td>2016</td><td>8</td><td>0</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Schleiz</td><td>Germany</td><td>50,57</td><td>11,81</td><td>2016</td><td>33</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Ziegelbach</td><td>Germany</td><td>47,89</td><td>9,86</td><td>2016</td><td>31</td><td>0</td><td>0</td><td>11</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Abercorn</td><td>Great Britain</td><td>56,00</td><td>-3,46</td><td>2014</td><td>11</td><td>0</td><td>0</td><td>6</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Abercorn</td><td>Great Britain</td><td>51,18</td><td>-3,47</td><td>2014</td><td>5</td><td>0</td><td>0</td><td>10</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Abercorn</td><td>Great Britain</td><td>56,00</td><td>-3,48</td><td>2014</td><td>1</td><td>0</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Abercorn</td><td>Great Britain</td><td>51,18</td><td>-3,47</td><td>2014</td><td>4</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Abercorn</td><td>Great Britain</td><td>56,00</td><td>-3,47</td><td>2014</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Abercorn</td><td>Great Britain</td><td>56,00</td><td>3,47</td><td>2014</td><td>2</td><td>0</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Blackmoor</td><td>Great Britain</td><td>51,18</td><td>-3,94</td><td>2015</td><td>30</td><td>1</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>London</td><td>Great Britain</td><td>51,51</td><td>-0,13</td><td>2010</td><td>27</td><td>4</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Thorpe</td><td>Great Britain</td><td>53,04</td><td>-1,78</td><td>2015</td><td>17</td><td>0</td><td>0</td><td>15</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Varybobi</td><td>Greece</td><td>38,13</td><td>23,78</td><td>2016</td><td>50</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Vestino</td><td>Greece</td><td>39,38</td><td>22,75</td><td>2016</td><td>112</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Brunico</td><td>Italy</td><td>46,80</td><td>11,96</td><td>2016</td><td>46</td><td>5</td><td>2</td><td>3</td><td>4</td><td>0</td><td>0</td></tr><tr><td>Mammola</td><td>Italy</td><td>38,38</td><td>16,19</td><td>2016</td><td>19</td><td>3</td><td>0</td><td>20</td><td>0</td><td>5</td><td>0</td></tr><tr><td>Teramo</td><td>Italy</td><td>42,64</td><td>13,67</td><td>2016</td><td>14</td><td>0</td><td>0</td><td>41</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Kalnciemi</td><td>Latvia</td><td>57,66</td><td>25,44</td><td>2016</td><td>165</td><td>0</td><td>26</td><td>46</td><td>0</td><td>0</td><td>3</td></tr><tr><td>Mazlimbas</td><td>Latvia</td><td>57,71</td><td>25,40</td><td>2016</td><td>75</td><td>0</td><td>16</td><td>65</td><td>2</td><td>4</td><td>0</td></tr><tr><td>Oubogo</td><td>Macedonia</td><td>41,63</td><td>21,99</td><td>2016</td><td>34</td><td>1</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Ait Oumribet</td><td>Morocco</td><td>29,54</td><td>-9,59</td><td>2016</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Bab Taza</td><td>Morocco</td><td>31,26</td><td>-5,19</td><td>2016</td><td>1</td><td>2</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Bounaamane</td><td>Morocco</td><td>29,54</td><td>-9,80</td><td>2016</td><td>0</td><td>43</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Chefchaouen</td><td>Morocco</td><td>35,16</td><td>-5,27</td><td>2016</td><td>18</td><td>4</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Dardara</td><td>Morocco</td><td>34,49</td><td>-5,30</td><td>2016</td><td>2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Moulay</td><td>Morocco</td><td>31,26</td><td>-6,25</td><td>2016</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Ouijane</td><td>Morocco</td><td>29,62</td><td>-9,54</td><td>2016</td><td>0</td><td>19</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Oulad Mtaa</td><td>Morocco</td><td>31,26</td><td>-8,26</td><td>2016</td><td>1</td><td>0</td><td>0</td><td>32</td><td>0</td><td>1</td><td>0</td></tr><tr><td>Ourtzagh</td><td>Morocco</td><td>34,49</td><td>-5,11</td><td>2016</td><td>0</td><td>3</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Rabat</td><td>Morocco</td><td>33,96</td><td>-6,80</td><td>2010</td><td>4</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Sidi Badhaj</td><td>Morocco</td><td>31,27</td><td>-8,22</td><td>2016</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Tizgzaouine</td><td>Morocco</td><td>31,77</td><td>-8,13</td><td>2016</td><td>1</td><td>3</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Tnine Aglou</td><td>Morocco</td><td>30,42</td><td>-9,76</td><td>2016</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Lomm</td><td>Netherlands</td><td>51,45</td><td>6,20</td><td>2017</td><td>18</td><td>0</td><td>0</td><td>10</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Olst</td><td>Netherlands</td><td>52,37</td><td>6,12</td><td>2017</td><td>11</td><td>0</td><td>0</td><td>6</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Narvik</td><td>Norway</td><td>68,48</td><td>17,51</td><td>2015</td><td>82</td><td>0</td><td>21</td><td>3</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Storforshei</td><td>Norway</td><td>66,31</td><td>14,81</td><td>2015</td><td>4</td><td>0</td><td>91</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>Tiller</td><td>Norway</td><td>63,34</td><td>10,42</td><td>2015</td><td>19</td><td>2</td><td>50</td><td>4</td><td>1</td><td>2</td><td>0</td></tr></tbody></table>
<table><tr><td>Jaksamanice</td><td>Poland</td><td>49,75</td><td>22,86</td><td>2016</td><td>93</td><td>2</td><td>22</td><td>26</td><td>0</td><td>3</td><td>0</td></tr><tr><td>Olszewnica</td><td>Poland</td><td>51,76</td><td>22,51</td><td>2016</td><td>25</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Pikulice</td><td>Poland</td><td>49,75</td><td>22,79</td><td>2016</td><td>195</td><td>1</td><td>23</td><td>29</td><td>0</td><td>2</td><td>0</td></tr><tr><td>Wronka</td><td>Poland</td><td>54,00</td><td>21,68</td><td>2016</td><td>0</td><td>0</td><td>35</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Aboadela</td><td>Portugal</td><td>41,27</td><td>-8,00</td><td>2009</td><td>26</td><td>2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Bajouca</td><td>Portugal</td><td>39,87</td><td>-8,77</td><td>2009</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Aleksandrovo</td><td>Serbia</td><td>45,63</td><td>20,58</td><td>2016</td><td>61</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Cirikovac</td><td>Serbia</td><td>44,67</td><td>21,19</td><td>2016</td><td>78</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Gamzigrad</td><td>Serbia</td><td>43,91</td><td>22,18</td><td>2016</td><td>71</td><td>1</td><td>0</td><td>3</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Pacir</td><td>Serbia</td><td>45,89</td><td>19,45</td><td>2016</td><td>59</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Arauut</td><td>Spain</td><td>41,97</td><td>1,76</td><td>2017</td><td>30</td><td>0</td><td>0</td><td>6</td><td>0</td><td>1</td><td>0</td></tr><tr><td>Arucas</td><td>Spain</td><td>28,13</td><td>-15,50</td><td>2010</td><td>37</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Castril</td><td>Spain</td><td>37,73</td><td>-2,34</td><td>2017</td><td>16</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Massanes</td><td>Spain</td><td>41,77</td><td>2,67</td><td>2017</td><td>0</td><td>0</td><td>0</td><td>17</td><td>0</td><td>20</td><td>0</td></tr><tr><td>Orgiva</td><td>Spain</td><td>36,85</td><td>-2,99</td><td>2017</td><td>6</td><td>0</td><td>0</td><td>4</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Palma</td><td>Spain</td><td>39,56</td><td>2,65</td><td>2015</td><td>26</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Torre del Campo</td><td>Spain</td><td>37,83</td><td>-3,90</td><td>2017</td><td>13</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Dittingen</td><td>Switzerland</td><td>47,44</td><td>7,50</td><td>2009</td><td>18</td><td>0</td><td>0</td><td>22</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Muhleturnen</td><td>Switzerland</td><td>46,81</td><td>7,51</td><td>2009</td><td>15</td><td>0</td><td>0</td><td>21</td><td>1</td><td>0</td><td>0</td></tr><tr><td>Kuzgöl</td><td>Turkey</td><td>40,91</td><td>32,10</td><td>2016</td><td>2</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Şahinköy</td><td>Turkey</td><td>40,24</td><td>28,24</td><td>2016</td><td>31</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></table>
<table><caption>**Tableau S2.** Primers used for PCRs and sequencing in this study.</caption><thead><tr><th>Target<br>Marker</th><th>Forward primer (5' to 3')</th><th>Reverse Primer (5' to 3')</th><th>Length of<br>amplified<br>fragment (bp)</th><th>References</th></tr></thead><tbody><tr><td>Cox1</td><td>LCO1490:<br>GGTCAACAAATCATAAAGATATTGG</td><td>HCO2198:<br>TAAACTTCAGGGTGACCAAAAAATCA</td><td>710</td><td>Folmer et al.<br>1994</td></tr><tr><td>16S<br>rDNA</td><td>16S-EK_F:<br>CACGTAAGAACTAAATAGTCGAAC</td><td>16S-EK_R:<br>GACCGTGCAAAGGTAGCATAATC</td><td>450</td><td>Ekrem et al.<br>2010</td></tr><tr><td>28S<br>rDNA</td><td>28S_S3660:<br>GAG AGT TMA ASA GTA CGT GAA<br>AC</td><td>28S_A335:<br>TCG GAA GGA ACC AGC TAC TA</td><td>657</td><td>Dowton &<br>Austin 1998<br>Whiting et al.<br>1997</td></tr></tbody></table>
<table><caption>Tableau S3. Reference sequences used for specific assignation.</caption><thead><tr><th>Species</th><th>Accession number</th><th>Origin</th></tr></thead><tbody><tr><td>Culicoides abchazicus</td><td>MK760081</td><td>GenBank</td></tr><tr><td>Culicoides abchazicus</td><td>MK760082</td><td>GenBank</td></tr><tr><td>Culicoides alachua</td><td>MK760088</td><td>GenBank</td></tr><tr><td>Culicoides chiopterus</td><td>KT278212</td><td>GenBank</td></tr><tr><td>Culicoides chiopterus</td><td>KM941505</td><td>GenBank</td></tr><tr><td>Culicoides chiopterus</td><td>KM940185</td><td>GenBank</td></tr><tr><td>Culicoides dewulfi</td><td>KJ162977</td><td>GenBank</td></tr><tr><td>Culicoides montanus</td><td>MK760206</td><td>GenBank</td></tr><tr><td>Culicoides montanus</td><td>MK760202</td><td>GenBank</td></tr><tr><td>Culicoides montanus</td><td>MK760207</td><td>GenBank</td></tr><tr><td>Culicoides montanus</td><td>MK760205</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus</td><td>DQ162811</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus</td><td>HM022839</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus</td><td>HM022842</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus</td><td>HM022856</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus</td><td>HQ824374</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus</td><td>KJ624104</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade dark</td><td>KR057928</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade dark</td><td>KR057935</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade dark</td><td>KR057940</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade dark</td><td>KR057929</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade dark</td><td>KR057930</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade dark</td><td>KR057932</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O1</td><td>HQ824371</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O1</td><td>HQ824378</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O1</td><td>JQ620132</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O2</td><td>HQ824381</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O2</td><td>HQ824382</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O2</td><td>HQ824383</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O2</td><td>HQ824384</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O3</td><td>JQ620144</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O3</td><td>JQ978454</td><td>GenBank</td></tr><tr><td>Culicoides obsoletus clade O3</td><td>JQ978459</td><td>GenBank</td></tr><tr><td>Culicoides sanguisuga</td><td>MK760237</td><td>GenBank</td></tr><tr><td>Culicoides sanguisuga</td><td>MK760238</td><td>GenBank</td></tr><tr><td>Culicoides scoticus</td><td>HM022865</td><td>GenBank</td></tr><tr><td>Culicoides scoticus</td><td>HM022875</td><td>GenBank</td></tr><tr><td>Culicoides scoticus</td><td>KJ624126</td><td>GenBank</td></tr><tr><td>Culicoides scoticus</td><td>KT186881</td><td>GenBank</td></tr><tr><td>Culicoides sinanoensis</td><td>MK760244</td><td>GenBank</td></tr><tr><td>Culicoides sinanoensis</td><td>MK760242</td><td>GenBank</td></tr><tr><td>Culicoides sp.</td><td>GMGRC1056-13</td><td>BOLD Systems</td></tr><tr><td>Culicoides sp.</td><td>GMGRC1000-13</td><td>BOLD Systems</td></tr><tr><td>Culicoides sp.</td><td>GMGRD2587-13</td><td>BOLD Systems</td></tr></tbody></table>
**Tableau S4.** Information on all *Culicoides* sequenced and GenBank accession numbers.

The table being too voluminous, it is not included in this manuscript.

DOI: https://doi.org/10.1186/s13071-020-04114-1
<figure><img src="image_31.png" /><figcaption>**Figure S1.** Maximum likelihood phylogenetic tree using 16S rDNA. Values at the nodes represent bootstrap values greater than 800 for the most important nodes (1000 replicates).</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>**Figure S2.** Maximum likelihood phylogenetic tree using *cox1*. Values at the nodes represent bootstrap values greater than 800 for the most important nodes (1000 replicates).</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>**Figure S3** -Maximum likelihood phylogenetic tree using 28S rDNA. Values at the nodes represent bootstrap values greater than 800 for the most important nodes (1000 replicates).</figcaption></figure>
Chapitre II : Recent phylogeographic history in
western Europe of one main species vector of
bluetongue and Schmallenberg viruses, *Culicoides*
*obsoletus* (Diptera: Ceratopogonidae) demonstrates
important continental connectivity
*Culicoides obsoletus* est une espèce autochtone et très abondante en Europe de l'ouest, mais la cartographie de la diversité cryptique a mis en évidence des variations de son abondance relative. *Culicoides obsoletus* est quasi omniprésent d dans la région des Balkans. Au contraire, on le trouve en plus faible proportion dans le nord de l'Europe et beaucoup plus rarement en Afrique du nord. Les épisodes glaciaires du Pléistocène ont modifié l'aire de répartition géographique, la structure génétique et l'histoire démographique de la plupart des espèces européennes. Si l'histoire évolutive de *C. imicola* a été décrite, elle reste la seule espèce pour laquelle des hypothèses approfondies ont été faites pour expliquer la distribution de sa diversité actuelle.

Ce deuxième chapitre présente une étude multi-marqueurs de génétique des populations et de phylogéographie (11 marqueurs microsatellites nouvellement développés et le gène mitochondrial Cox1) menée à l'échelle européenne sur *C. obsoletus* pour identifier la structure génétique et phylogéographique de l'espèce et déterminer les événements historiques qui ont façonné le schéma actuel de la structure génétique de cette espèce.

Nos travaux permettent pour la première fois de déterminer, en observant les schémas de distribution des fréquences alléliques, comment la dispersion et les événements historiques ont pu façonner l'agencement spatial de la diversité génétique. Notre étude a révélé de faibles niveaux de différenciation génétique entre les populations à l'échelle européenne. Les méthodes bayésiennes de classification révèlent un cluster principalement représenté dans la zone des Balkans, un représentatif des îles Canaries, et un autre englobant le reste des populations Européennes. Nous avons montré une absence de signature phylogénétique et l'analyse du réseau d'haplotypes met également en évidence un manque de structuration spatiale. Nos résultats soutiennent une expansion démographique récente. Les méthodes bayésiennes ABC nous ont permis de mettre en évidence une divergence récente (au cours des 10 000 dernières années) des populations de *Culicoides*, avec une recolonisation post-glaciaire à partir des refuges nordiques et balkaniques, suivie d'une expansion géographique vers l'ouest de l'Europe des populations issues des refuges, créant des populations mélangées à l'ouest de l'Europe.

Ce chapitre présente des résultats préliminaires qui sont encore à affiner avant soumission pour publication dans *Scientific Reports*.
Recent phylogeographic history in western Europe of one main species vector
of bluetongue and Schmallenberg viruses, *Culicoides obsoletus* (Diptera:
Ceratopogonidae) demonstrates important continental connectivity

Antoine Mignotte<fnref n="1" />,<fnref n="2" />*, Claire Garros<fnref n="1" />,<fnref n="2" />, Laetitia Gardès<fnref n="1" />,<fnref n="3" />, Thomas Balenghien<fnref n="1" />,<fnref n="2" />,<fnref n="4" />, Maxime Duhayon<fnref n="1" />,<fnref n="2" />, Ignace Rakotoarivony<fnref n="1" />,<fnref n="2" />, Karine Huber<fnref n="1" />.

<footnote n="1">¹ ASTRE, Univ Montpellier, Cirad, INRA, Montpellier, France.</footnote>

<footnote n="2">² Cirad, UMR ASTRE, F-34398 Montpellier, France.</footnote>

<footnote n="3">³ Cirad, UMR ASTRE, F-97170 Petit-Bourg, Guadeloupe, France.</footnote>

<footnote n="4">⁴ Institut Agronomique et Vétérinaire Hassan II, Unité Parasitologie et maladies parasitaires, 10100 Rabat, Morocco.</footnote>

<footnote>* Corresponding author: Antoine Mignotte, UMR 117 ASTRE, Campus international de Baillarguet, 34398 Montpellier cedex 5, antoine.mignotte@cirad.fr; or Claire Garros, UMR 117 ASTRE, Campus international de Baillarguet, 34398 Montpellier cedex 5, claire.garros@cirad.fr</footnote>
Since 2006, Northern Europe faced massive outbreaks of bluetongue virus (BTV) which affects
wild and domestic ruminants. This was followed by the emergence of Schmallenberg virus
(SBV) in 2011, another Culicoides-borne virus, which also caused important economic losses
to livestock breeders. Culicoides obsoletus has been identified as the main vector species in
western Europe. Culicoides obsoletus is autochtonous and very abundant in Europe, in
particular in the Balkan Peninsula. On the contrary, it is found at a lower proportion in northern
Europe and much rarer in North Africa. During the Pleistocene glacial episodes, the
geographical range, genetic structure and demographic history of most European species was
impacted. Moreover, Culicoides dispersion is referred to as "stratified dispersion", which results
from the combination of processes occurring at short distances and jumps over long distances.
The objective of this work is to determine, by observing patterns of allelic frequency
distributions, how dispersion and historic events may have shaped the spatial arrangement of
genetic diversity. This study is conducted at continental scale, on two transects from Morocco
to Finland and from Turkey to the United Kingdom. We used 11 microsatellite markers and a
mitochondrial gene (cox1), to identified genetic and phylogeographic structure of C. obsoletus
at the European scale and determine historical events that have shaped the current genetic
structure pattern. Our study highlights important gene flows at the continental level and the
connectivity between populations separated by a high geographical distance. The use of
mitochondrial and microsatellite markers has allowed us to highlight a lack of phylogeographic
and genetic structure, a recent increase in effective population size link with post-glacial
recolonization from Nordic and Balkanic refuges, then a geographical expansion towards the
West of Europe of the populations coming from the refuges, creating populations admixed in
the West of Europe. The recent population expansion highlights the potential role of the glacial
episode and the post-glacial colonization event in shaping the geographic distribution of C.
obsoletus.

Key words: *Culicoides obsoletus*, ABC, Dispersion, Microsatellites, Genetic structure,
Palearctic region, population genetics, Phylogeography
Introduction

In 2006, western Palearctic region faced outbreaks of bluetongue virus (BTV) followed by Schmallenberg virus (SBV) in 2011 (Carpenter et al, 2013; Carpenter et al, 2009; Hoffmann et al, 2009; Savini et al, 2008). The unexpected nature of 2006 outbreaks compared to the outbreaks in the Mediterannean basin in the 1990-2000's is based on the absence of C. imicola in the western Palearctic region (Carpenter et al, 2009; Mellor et al, 2008; Purse et al, 2005). Massive spread in the western Palearctic region raises the question on the vector competence of native Culicoides (Carpenter et al, 2009; Wilson and Mellor, 2009). Among them, Culicoides obsoletus has been highlighted as the main vector species (Carpenter et al, 2009; Hoffmann et al, 2009; Savini et al, 2008).

*Culicoides obsoletus* is an native and abundant species in Europe and able to breed in a wide range of habitats (Mehlhorn *et al*, 2007; Purse *et al*, 2015). Initially described in Europe, many synonymous species have been described as *varius* ((Winnertz), 1852) in Germany, *yezoensis* ((Matsumura), 1911) in Russia, *lacteinervis* (Kieffer, 1919) in Ukraine, *kabyliensis* (Kieffer, 1922) in Algeria, or *intermedius* (Okada, 1941) in Japan. *Culicoides obsoletus* has a wide latitudinal distribution from North Africa (Morocco (Bourquia et al, 2019), Tunisia (Sghaier et al, 2017), Algeria (Hakima et al, 2020)) to northern Europe (Finland and Sweden (Mignotte et al, 2020)). It is also found on islands such as the Balearic Islands (Miranda et al, 2003), the Canary Islands (Martinez-de la Puente et al, 2012), and Corsica (Delecolle, 1985). It is thus mainly present in the northern hemisphere (Borkent, 2016). However, although its distribution center seems to be in Europe, it is also present on two other areas of the northern hemisphere (Borkent, 2016), America (Canada (Anderson et al, 1991) and USA (Borkent, 2016), and Asia (China (Gao et al, 2020), South Korea (Cho and Chong, 1974)). This species is therefore considered to be Palearctic and Nearctic (Borkent, 2016; Mathieu, 2011). Its trophic behaviour is flexible, allowing it to be present in a wide range of habitats. In addition, *C. obsoletus* is very generalist and is capable of blood feeding on a wide range of hosts such as goats, horses, sheep, rodents and birds (Garros and Balenghien, 2017). The abundance of *C. obsoletus* varies greatly depending on the country (Cuellar *et al*, 2018). A recent study describing the diversity within the Obsoletus/Scoticus Complex showed a pattern of decreasing proportion of *C. obsoletus* numbers captured following an East-West European gradient (Mignotte *et al*, 2020). Indeed, in the Balkan Peninsula, *C. obsoletus* is the most abundant species of the complex almost the only species caught during the sampling surveys. This has also been highlighted in Romania,
Slovakia, Albania, Macedonia, Kosovo, Montenegro, Serbia, Bosnia, Bulgaria and Croatia (Goffredo *et al*, 2004; Hristescu *et al*, 2020; Mignotte *et al*, 2020; Pudar *et al*, 2018; Sarvasova *et al*, 2014). On the contrary, it is found at a lower proportion in northern Europe (Mignotte *et al*, 2020) and much rarer in North Africa (Bourquia *et al*, 2019; Hakima *et al*, 2020; Sghaier *et al*, 2017). In addition, *C. obsoletus* has a strong intra- and interspecific diversity distributed throughout Europe, with western areas, such as France and Italy, where all the members of the Obsoletus/Scoticus Complex have recently been identified (Mignotte *et al*, 2020).

Geographic range, genetic structure and demographic history of most European species was modified during the glacial episodes of the Pleistocene (Bowler and Benton, 2005; Hewitt, 2004; Phillips et al, 2010; Widmer and Lexer, 2001). As a result of the rise in temperatures during interglacial episodes, new ecological niches became vacant in ice free regions, allowing the recolonization of no cold adapted species from glacial refuges (Hewitt, 2004; Petit et al, 2003; Taberlet et al, 1998). This post-glacial colonization also allows the creation of reconnection areas. In Europe, the Balkan, Italian and Iberian peninsulas have played the role of glacial refuges due to the reduction of viable areas during the glacial episodes (Hewitt, 1999; Petit et al, 2003; Taberlet et al, 1998). During the interglacial periods, these peninsula were the geographical origin of the spatial and demographic expansion of many species (Abellan and Svenning, 2014; Gómez and Lunt, 2007; Rodríguez-Sánchez et al, 2010). For example, the Balkan Peninsula has served as glacial refuges for butterflies (Schmitt et al. 2006a). In addition, recent studies show that in some areas in central and eastern Europe with stable microclimates may also have acted as glacial refuges for some cold tolerant species (Bhagwat and Willis, 2008; Schmitt and Varga, 2012; Sommer and Zachos, 2009; Stewart and Lister, 2001; Whittaker and Field, 2000).

The dispersal capacity of insect vectors is an important factor in the development of epidemiological models for predicting the evolution of epidemics. Indeed, dispersal capacity of insect vectors can be used to identify spatial pathogen transmission pattern (Kozakiewicz et al, 2018). The dispersion of *Culicoides* corresponds to what ecologists call 'dispersive stratified' which results from the combination of processes occurring at different spatial scales, that is to say combination of dispersal processes occurring actively at short distances, and passively or semi actively with jumps over long distances with wind-mediated dispersal patterns (Jacquet et al, 2016a; Murray and Kirkland, 1995; Onyango et al, 2015a). Important gene flows have been demonstrated between remote geographical areas in *Culicoides* species in Mediterranean basin
(Jacquet *et al*, 2016a), Africa (Onyango *et al*, 2015a) and in Australia (Pili *et al*, 2010). Dispersal capacity is an important component during periods of recolonization of newly available habitats, as is the case during inter-glacial periods. This allows for rapid spatial and demographic expansion towards ecological niches that have become vacant again.

Gene cytochrome oxidase I (cox1) has been used in phylogeographic studies on *C. imicola* in order to identify the colonization routes of this initially Afrotropical vector in the Mediterranean basin (Calvo *et al*, 2009; Nolan *et al*, 2008). Maternally inherited, the mitochondrial DNA has a relatively high rate of evolution and does not recombine between loci, which makes it a coherent marker for reconstructing the pathways of taxa during their recolonisation (Estoup and Guillemaud, 2010; Kirk *et al*, 2013). Microsatellite markers are complementary for studying population connectivity and demographic events (Hewitt, 2004). The use of a multi-local approach allows for a more refined and time-scaled estimate of connectivity between populations (Brito and Edwards, 2009; Prugnolle and de Meeus, 2002). We characterized the genetic structure using Bayesian clustering and traditional population genetics tools Approximate Bayesian Computation (ABC) (Beaumont and Rannala, 2004) were used to identify evolutionary scenarios to retrace and date demographic and historical events shaping the current population genetic structure of *C. obsoletus*. The objective of our study is to identify species connectivity on a large geographical scale of main species vector of bluetongue and Schmallenberg viruses in Europe. For this, we used 11 microsatellite markers and a mitochondrial gene (cox1), to identify genetic and phylogeographic structure of *C. obsoletus* at the European scale and determine historical events that have shaped the current genetic structure pattern. We used populations from 20 countries along a Palaearctic Mediterranean transect covering Scandinavia to Canary Islands (North to South) and Canary Islands to Turkey (West to East). These results will allow us to clarify the study of the dispersal of *C. obsoletus* on land, but also to know the phylogeographic history of one main vectors of livestock virus in Europe.
# Methods

## Sampling and species morphological identification

We collected a total of 2409 individual female biting midges in 19 countries, from 59 sampling sites in the Palearctic region between 2009 and 2017 using national surveillance networks for *Culicoides* populations or local collections (Table S1). Collection and morphological identification were carried out according to the same protocol as in a previous study on the phylogeny of *Culicoides* (Mignotte *et al*, 2020).

## DNA extraction, amplification and sequencing

A total of 2409 individuals belonging to the *C. obsoletus* species complex were DNA-extracted, cox1 amplified and sequenced following the procedure detailed in Mignotte *et al.* (2020). To ensure that only *C. obsoletus* individuals were included in the final data set, a DNA barcode analysis was performed according to Mignotte *et al.* (2020). All cox1 sequences were then aligned using the MUSCLE algorithm (Edgar, 2004) implemented in the software package GENEIOUS v.6.0.5 (Biomatters, http://www.geneious.com). In order to estimate unbiased gene flows, only individuals identified as *C. obsoletus* were used for further analyses.

Among these individuals, 829 were amplified and genotypes (Tab S6) according to (Mignotte *et al*, Submitted in 2020). Allelic size allocations for all individuals and microsatellite markers were performed using the program GeneMapper® version 5 (Life Technologies, Applied Biosystem) with double blind reading to limit the potential interpretation bias from the reader.

## Mitochondrial analyses

We calculated nucleotide diversity ($\pi$), haplotype number (H) and diversity (Hd) using the program DNASP v.5.10 (Librado and Rozas, 2009). In order to detect a past demographic change we performed neutrality tests based on Tajima's D and Fu's Fs statistics with DNASP v.5.10. (Fu, 1997; Ramos-Onsins and Rozas, 2002; Tajima, 1989). In order to detect sites under adaptive selection, we run fixed-effects likelihood (FEL) and single-likelihood ancestor counting (SLAC) on the cox1 data set. These analyses were run in the DataMonkey server (http://www.datamonkey.org) (Kosakovsky Pond and Frost, 2005; Murrell *et al*, 2012). We inferred demographic history using mismatch distribution analysis (Rogers and Harpending, 1992) with ARLEQUIN 3.5 program. Due to the lack of spatial phylogeographic structuring, it was not possible to define a similar population grouping so all populations were therefore treated as a single population.
## Genetic structure, diversity and dissimilarity

Matrix of pairwise $\Phi_{ST}$ among populations (Excoffier *et al*, 1992) was calculated with SPADS, and used to investigate mitochondrial isolation-by-distance. In order to visualize genealogical relationships among European populations, a median-joining haplotype network of cox1 mitochondrial sequences of *C. obsoletus* was constructed with POPART (REF). Mitochondrial data set was used to evaluate population differentiation $G_{ST}$ (Pons and Petit, 1995) that makes use only of the allelic frequencies, $N_{ST}$ (Pons and Petit, 1996) for which similarities between the haplotypes are considered. The phylogeographic structure is evaluated by testing $N_{ST} - G_{ST}$ (Pons and Petit, 1996) with 10,000 permutations using SPADS (Dellicour and Mardulyn, 2014). Program MAPI was used based on 1,000 permutations and an $\alpha$ value of 0.05, to map inter-population genetic distances, $\Phi_{ST}$, to visualize the spatial variation in genetic dissimilarity and define areas where genetic dissimilarity is significantly lower or higher than expected by chance (Piry *et al*, 2016).

## Microsatellites analyses

### Population structure and genetic diversity

We calculated number of alleles per locus, the allele size range, the observed and expected heterozygosity were estimated per locus, and per population over the entire dataset using the R packages *poppr* (Kamvar *et al*, 2014) and *mmod* (Winter, 2012). We calculated in R the Spearman's rank correlation coefficient between the allele sizes and their frequencies (Wattier *et al*, 1998) to identify a potential bias related to short allele dominance (Manangwa *et al*, 2019; Wattier *et al*, 1998). Null allele frequencies were calculated for each pair of loci using the methods of Chakraborty *et al* (1994) and Brookfield (1996) (Brookfield, 1996; Chakraborty *et al*, 1994) implemented in R package *PopGenReport* (Adamack *et al*, 2014). Linkage disequilibrium between each pair of loci were performed by the calculation of index of association (Agapow and Burt, 2001) implemented in R package *poppr* (Kamvar *et al*, 2014). Bonferroni's correction was applied to linkage disequilibrium analyses to take into account potential biases associated with multiple comparisons (Rice, 1989). We calculated the fixation index ($F_{IS}$) and tested the deviation from the Hardy-Weinberg equilibrium, (Weir and Cockerham, 1984) using FSTAT v2.9.3.2 (Goudet *et al*, 1996). Correlation between $F_{IS}$ and null allele frequencies was calculated to try to explain a potential deviation from the Hardy-Weinberg equilibrium.
We used the program STRUCTURE (version 2.3.4) (Pritchard *et al*, 2000) to estimate the most probable number of genetic clusters (K) and individual probability of assignment to each genetic cluster. In order to assume identical allele frequencies among populations, Bayesian clustering analysis was run with correlated allele frequency. Admixture model was used to estimate for each individual the belonging proportion to each K clusters. Particularly adapted when there is a weak genetic structuring (Hubisz *et al*, 2009), the “locprior” model was used to use sampling locations as prior information to assist the clustering. Value of tested K varying from 1 to 10, with 10 independent runs, with a burn-in of 100,000 iterations in the Markov Monte Carlo chain (MCMC) and a stationary phase of 1,000,000 iterations in the MCMC. We used the STRUCTURE Harvester web platform (Earl and vonHoldt, 2011) to identify most probable number of clusters with the ΔK method (Evanno *et al*, 2005). In order to visualize the partition in genetic differentiation among all samples a neighbor-joining (NJ) tree (Takezaki and Nei, 1996) was constructed with the software POPULATIONS v.1.2.30 (http://bioinformatics.org/~tryphon/populations/) from a Cavalli-Sforza and Edwards genetic distance (1967) matrix. The trees were then subjected to 1000 bootstrap replications to confirm topological accuracy. In order to investigate isolation-by-distance, we used the R package *fields* (Nychka *et al*, 2017) to perform Mantel tests with the R package *vegan* (Oksanen *et al*, 2013) and linear regressions between inter-population ($F_{ST}/(1-F_{ST})$) genetic distances and log-transformed geographic distances between sampling sites with R package *fields* (Nychka *et al*, 2017). Allelic richness, classically defined as the number of alleles present in a population at a given locus, was calculated and corrected for variations in sample sizes between populations using the rarefaction method of the R package *Hierfstat* (Weir and Goudet, 2017; Yang, 1998). We used NE ESTIMATOR software and linkage disequilibrium method to evaluate the effective population size of *Culicoides* populations using the microsatellite dataset (Do *et al*, 2014). The program MAPI was used, based on 1,000 permutations and an $\alpha$ value of 0.05, to map inter-population genetic distances (pairwise $F_{ST}$) to visualize the spatial variation in the genetic dissimilarity and to define areas where the genetic dissimilarity is significantly lower or higher than expected by chance (Piry *et al*, 2016).

## Scenarios inference with Approximate Bayesian computation

In order to infer evolutionary scenarios of *C. obsoletus* range expansion in Europe, we used Approximate Bayesian computation (ABC) (Beaumont *et al*, 2002; Bertorelle *et al*, 2010) approach using DIYABC software v.2.0.4 (Cornuet *et al*, 2014). ABC analyses were run on our microsatellite dataset. We used the population independently without grouping them together
as representative of the geographical area. The list of the mutational, demographic and historical prior of model parameters set implemented in DIYABC are present in table S3. We used mean number of alleles, mean allelic size variance, mean expected heterozygosity (Nei *et al*, 1975), pairwise FST values (Weir and Cockerham, 1984), mean ratio of the number of alleles over the range of allele size (Garza-Williamson's M) (Garza and Williamson, 2001), and the classification index (Pascual *et al*, 2007) as summary statistics. We assumed 4 generations per year for *C. obsoletus* (Mathieu, 2011). Microsatellite mutation rate use for ABC analysis ranging from $6 \cdot 10^{-6}$ to $10^{-4}$. Each scenario was tested using 1 million of simulated data sets. In order to identify the most probable scenario, we calculated the posterior probabilities associated with each scenario using a polychotomous logistic regression (Cornuet *et al*, 2014). This regression was performed on the 1% of the simulated data sets closest to the observed data set.

In order to reconstruct evolutionary history of *C. obsoletus* in Europe, different ABC scenarios were tested. This set of scenarios is aimed at identifying the pattern of divergence and the genealogy between the populations. For this purpose, the populations of Velestino (Greece), Narvik (Norway), Les Hitres (France) and Torre del Campo (Spain) were used. The populations included in this analysis are some of the most genetically differentiated, populations and are located in distant geographical areas, making it possible to trace phylogeographical history of *C. obsoletus* in Europe and test the hypothesis of European glacial refuges. All the tested scenarios are presented in table S2. In view of the main European glacial refuges and the genetic structure observed in the Mediterranean basin, the ABC scenarios were constructed to identify the geographical origin of post-glacial expansion at different dates of divergence (old: 20,000-100,000 years, recent: 5,000-25,000 years, or very recent: 2500 - 10,000 years). Four coalescence-based population evolution hypotheses on the demographic and phylogeographic history of *C. obsoletus* in Europe were tested with second set of scenario and different divergence time: (i) post-glacial recolonization from Nordic and Balkanic refuges, then a geographical expansion towards the West of Europe of the populations coming from the refuges, creating populations admixed in the West of Europe (ii) post-glacial recolonization from Iberic and Balkanic refuges, then a geographical expansion towards the North of Europe of the populations coming from the refuges, creating populations admixed in the North of Europe (iii) an origin from the Iberian Peninsula with a settlement from West to East (iiii) an origin from the Balkan Peninsula with a settlement from East to West (Figure 1).
# Results

## Mitochondrial analyses

We sequenced a fragment of 627 bp of the cox1 gene from 2409 *C. obsoletus* from 59 sampling sites (Table S1). We identified 106 different haplotypes in the cox1 dataset. The haplotype diversity ($H_d$) was 0,772 for all mitochondrial data set and a total of 67 variable sites out of 627 bp was observed. Haplotype diversity ranging from 0,895 (Dittingen, Switzerland) to 0,107 (Arucas, Spain). All the metrics showing the genetic diversity for each population of the cox1 dataset are shown in Table S4.

## Genealogical relationships, genetic structure and genetic diversity

The median-joining network based on cox1 sequences (Figure 2) shows, four widespread dominant haplotypes (H26, H86, H44, H5) grouping together individuals collected in distant geographical areas at the European scale. They are surrounded by several rare haplotypes distant from a few mutations. The median-joining network shows a star-like pattern signature characteristic of populations that have undergone a demographic expansion. The AMOVA $\Phi_{ST}$ estimator shows a significant genetic differentiation over all sampling sites ($\Phi_{ST} = 0.083$, p-value < 0.05). No significant phylogeographic signal was observed ($N_{ST} - G_{ST} = -1.357$, p-value = 1.0).

The mitochondrial genetic dissimilarity results from MAPI analyses show a relatively homogeneous pattern of dissimilarity, except for populations from Canary and Morocco regions. The maximum inter-individual genetic dissimilarity is of 0.495 and the lowest is of -0.240. The Canary Islands is the most genetically dissimilar region. On the contrary, the continental Europe together with British islands showed a strong genetic homogeneity (Figure 3.a). A significant but low isolation-by-distance was detected with the population genetic distance $\Phi_{ST}$ for the mitochondrial dataset ($r=0.020$, p=0.002, Figure 3.b). Mitochondrial nucleotide diversity (Pi) shows a heterogeneous distribution of diversity (Fig S3).

## Demographic history of *C. obsoletus* in western Palearctic Europe

The mismatch distributions plot showed a unimodal shape indicating recent demographic expansions in Europe (Rogers and Harpending, 1992) (Figure 4). We rejected neutrality for the neutrality analyses on the whole cox1 dataset with significant negative values observed for
Tajima's D (D = -2.17897, p-value < 0.05) and Fu's (D = -25.81725, p-value < 0.05) indicating a strong demographic expansion or natural positive selection events. No sites under positive selection were detected by SLAC and FEL analyses, and only 22 sites of 627 sites of cox1 were under negative selection (P value < 0.1). The effective population sizes were estimated for each population with Ne Estimator and range from 17.1 to infinite value (Table S7).

## Microsatellites analyses

*Microsatellites clustering analyses, mapping of genetic distances and genetic diversity*

We genotyped 829 *C. obsoletus* from 27 sampling sites located in the western Palearctic region (Table S6) at 11 microsatellite loci (Table S5). The number of alleles per locus varied from 14 (OBSms29) to 47 (OBSms5) for the whole dataset (Supplementary Table 8). All populations are at Hardy-Weinberg equilibrium (Supplementary Table 7). All population genetic diversity metrics are shown in the table S7. FIS values per locus are significantly correlated with null allele values (r=0.953, p-value>0.001). No significant correlation between allele size and allele frequency was found for any of the 11 loci (Table S7). The index of association does not revealed any linkage disequilibrium between all loci after Bonferroni correction for multiple comparisons.

The Bayesian clustering analysis performed with STRUCTURE on microsatellite data identified four genetic groups (ΔK maximum for K=4) (Figure S1). Cluster 1 is the majority in almost all the populations of Continental Europe and the British Isles (Figure 5). The populations of Continental Europe revealed a low genetic structure of *C. obsoletus*. The second cluster is the majority in two populations of Greece and Turkeys. The third cluster is only represented in majority in the Canary Islands and in Serbia. The results for the second optimal K value (K=2) of the STRUCTURE program gave a similar result, reinforcing the conclusion that the populations are genetically weakly structured, with only two populations from Greece and Turkey standing out (Figure S2). The rarefied allelic richness of the microsatellite dataset varies from 5,903 (Canary Islands) to 8,513 alleles (Turkey). The Balkan region is the region with the greatest allelic richness (Macedonia and Turkey). On the contrary, the populations of Northern Europe and Western Europe are the ones with the lowest allelic richness. No clear geographical structure was identified between all sample using neighbor-joining tree based on Cavalli-Sforza and Edwards genetic distance (FigS4).
The microsatellite genetic dissimilarity results from MAPI analyses showed a non-homogeneous pattern of dissimilarity. The maximum inter-individual genetic dissimilarity was of 0.042 and the lowest was of 0.011. The southern Europe was the most genetically dissimilar region. On the contrary, northern Europe showed low genetic dissimilarity (Figure 6.a). However, almost none of the areas was significantly genetically dissimilar in MAPI analysis with the exception of a few areas in northern Europe showing significant genetic dissimilarity. A significant isolation-by-distance was detected with the population genetic distance FST/(1-FST) for microsatellite dataset (R²=0.084, p=0.003, Figure 6.b).

## Scenarios inference with Approximate Bayesian Computation

The posterior probabilities calculated for each scenario provided higher statistical support to the scenarios 1 (p = 0.6350) (Table 1). This scenario assumes an origin from the Balkan peninsula and from Eastern Europe followed by an admixture event, with a recent divergence time (10,000-2500 years) (Table S2). All the tested scenarios and results are presented in table S2. Type I and II errors were calculated in order to evaluate to what extent this scenario could be trusted. Type I error rate (proportion of simulations in which the scenario considered as true is not true) for true scenario was 0.432. Type II errors (proportion of simulations in which the true scenario is not considered true) are 0.142 for scenarios 1. These relatively high values of error could reflect that our data are not enough informative or that our scenario could be more complex. The real and simulated data are relatively well grouped by Principal Component analysis (PCA) (Fig. S6), indicating that some priorities could be refined, but that our scenario is still consistent in explaining the real dataset. We can conclude that scenario 1 is the most likely scenario among the proposed scenarios provide relatively good description of observed dataset, but it will have to be completed.
Discussion

Combined approaches of population genetics and phylogeography, based on microsatellite and
mitochondrial markers have allowed us to explore the genetic structure but also the evolutionary
and demographic history of *C. obsoletus* at the European scale. Our study revealed low levels
of genetic differentiation between populations at the European scale. Bayesian clustering
methods reveal four genetic clusters distributed all over the study area but one of the clusters is
mainly represented in the Balkan area and another one encompassed the Canary Islands. We
have shown an absence of phylogenetic signature using the difference between N<sub>ST</sub> and G<sub>ST</sub>.
The haplotype network analysis highlights also a lack of spatial structuring. The significant
negative values of Fu's and Tajima's D support a strong demographic expansion. Approximate
Bayesian computation methods allowed us to highlight a recent divergence (10,000-2,500
years) of *Culicoides* populations, with a post-glacial recolonization from Nordic and Balkanic
refuges. Then, a geographical expansion towards the West of Europe of the populations coming
from the refuges, creating populations admixed in the West of Europe.

***A strong dispersal capacity that blurs the genetic structure***

The study of mitochondrial polymorphism does not reveal any genetic structure in our study area. Microsatellite markers do not make it possible to highlight a very marked genetic structure although a few slightly differentiated Mediterranean populations are observed. Mitochondrial-nuclear discordances are regularly reported in studies using a multi-markers approach (Kolleck et al, 2013; Rasic et al, 2015). This type of discordance is expected because the mitochondrial genome is haploid and uniparentally inherited in most animal species, and therefore has a fourfold smaller effective population size. Except incomplete lineage sorting, mitonuclear discordance can also arise if there are differences in how selection acts on the mitochondrial genome as compared to the nuclear genome or if there is biased movement of either marker type driven by demographic asymmetries, such as sex-biased dispersal (Frank et al, 2011). If selection for mtDNA variants varies geographically, then discordant patterns between mtDNA and nuDNA can arise (Irwin et al, 2005). Our results do not show a positive selection signal, but only a few sites under negative selection on mtDNA. However, these negatively selected sites are found in all sampling areas, which rules out selection as an explanatory factor of the observed spatial genetic patterns (Murrell et al, 2012). Female-biased dispersal can also promote the dispersal of mtDNA in the absence of concordant movement of nuDNA (Funk and Omland, 2003). If the dispersion of *Culicoides* by the wind on long distance is independent of
the sex it is not the case of the active dispersion on shorter distances. Indeed, in blood-feeding insects, females generally travel longer distances than males in search of egg-laying sites and blood meals (Kirkeby *et al*, 2013). The combination of these different mechanisms could therefore explain the slight signal differences between the two types of markers.

Only the population of the Canary Islands is well genetically differentiated. This population also exhibits low allelic richness and a very low nucleotide diversity. Island species are likely to have experienced population bottlenecks at some point in their evolutionary history due to founder events during the initial island colonization. As only a fraction of individuals from the original population found an island population, only a fraction of the original genetic diversity of the population will be maintained, and effective population sizes (*Ne*) will be small (James *et al*, 2016; Nei *et al*, 1975). As a result, this can lead to a reduction in genetic diversity and rapid genetic differentiation between these species and their native populations.

In many species of flying insects, there are important dispersal phenomena marked by a weak genetic structure (Lecocq *et al*, 2013; Lecocq *et al*, 2016; Mayer *et al*, 2015). For example, in Hymenoptera species of the genus *Bombus* a genetic differentiation is only observed between continental and island populations whereas no significant differentiation was found among populations from the European continent (Ellis *et al*, 2006; Estoup *et al*, 1996; Goulson *et al*, 2011; Lecocq *et al*, 2013; Shao *et al*, 2004; Widmer *et al*, 1998). This high level of gene flow is consistent with what we have observed for *C. obsoletus* at the scale of France (Mignotte *et al*, Submitted in 2020). This has also been reported for other *Culicoides* species in other geographical areas, in Australia for *C. brevittarsis* (Onyango *et al*, 2015a), in America for *C. stellifer* (Shults *et al*, 2020) and in Africa for *C. imicola* (Jacquet *et al*, 2016a).

We observed a significant but low level of isolation by distance in continental Europe that suggests a stepping-stone model of migration compatible with the combination of passive and active dispersal capacity of *Culicoides* species as it has already been demonstrated for *C. imicola* (Jacquet *et al*, 2016b). The passive dispersion of the *Culicoides* by winds over long distances certainly contributes to blurring the signal of isolation by distance. The intensification of grazing practices could allow to increase the size of *Culicoides* populations and facilitated the spread of BTV. The combined effects of active dispersion of *Culicoides* at short distances, passive dispersion by winds allowing long-distance jumping, and the recent intensification of breeding activities has certainly contributed to genetic homogenization, masking the past phylogeographic structure. *Culicoides obsoletus* is a species that feeds on a wide range of hosts and can lay eggs in a wide range of breeding sites (Garros and Balenghien, 2017). The generalist
and opportunistic behavior of this species, allowing them to adapt to a wide range of habitats
(Harrup *et al*, 2016), has also contributed to this genetic homogenization.

*A glacial refuge in the Balkan Peninsula during the Pleistocene*

The low species diversity within the Obsoletus/Scoticus Complex observed in southeastern Europe and the Balkans (Mignotte *et al*, 2020), combined with the high allelic richness observed for *C. obsoletus* in this region, led us to hypothesize that this area could have served as a refuge for this species during the Pleistocene. The Balkan peninsula have played the role of glacial refuge for many species due to the reduction of viable areas during the glacial episodes (Hewitt, 1999; Petit *et al*, 2003; Taberlet *et al*, 1998). The lack of gene flow between isolated populations due to the glaciation has allowed genetic differentiation of populations due to genetic drift and the selection process, leading to potential speciation events (Hewitt, 1999; Hewitt, 2004). The spatial expansion of populations from glacial refuge has been accompanied by a gradient of reduction in the genetic diversity of the population during the colonization process (Comes and Kadereit, 1998), caused by serial founder events, creating genetic bottlenecks (Dlugosch and Parker, 2008; White *et al*, 2014). The genetic diversity is on average higher in areas that have long functioned as refuges than in areas of post-glacial redeployment (Prugnolle *et al*, 2005; White *et al*, 2013). In the same way for species of eastern Palearctic origin, the existence of refuges located within the limits of the ice caps where many plant species survived the last glacial maximum, have been shown (Abbott *et al*, 2003; Lagercrantz and Ryman, 1990). The Approximate Bayesian computation methods allowed us to identify as the most probable scenario a scenario identifying the Balkans and the arctic region as the post-glacial origin of ancestral populations, 5,000 to 10,000 years ago, which after more recent admixture events (2,500 years ago) gave rise to the populations of Western Europe.

The DIYABC software used to conduct these analyses does not allow to include recurrent
migration phenomena in the tested scenarios (Cornuet *et al*, 2014), which can be a hindrance to
study the evolutionary history of species with high dispersal capacities such as *C. obsoletus*. A
high dispersal capacity can sometimes even prevent the identification of invasion routes, as was
the case during the rapid and accidental invasion in Europe of *Leptoglossus occidentalis*, an
insect native to North America (Lesieur *et al*, 2018). Moreover, some studies have suggested
the existence of refugia not only in southern Europe, but also in western, central and eastern
Europe (Schmitt and Varga, 2012; Stewart and Lister, 2001). *Culicoides obsoletus* has a very
wide ecological niche allowing it to adapt to a wide variety of environments. In consequence,
during episodes of glaciation, more scattered and northerly glacial refuges were able to allow
the persistence of coniferous and deciduous forests the animal species associated with them and
among it potentially *C. obsoletus* (Provan and Bennett, 2008; Tzedakis *et al*, 2013). *Culicoides obsoletus* can be found from sea level to high altitudes. This wide altitudinal distribution has
potentially allowed it, as it has already been shown for many species of boreal-mountain, a
broader distribution during glaciations than in interglacial periods, limiting the genetic and
phylogeographic structure of their populations (Schmitt and Varga, 2012; Taberlet *et al*, 1998).
Thus, to have a global vision of the phylogeographic history of the species it would be necessary
to integrate individuals from Asia and America. More locally, to dress a complete view of the
evolutionary history of *C. obsoletus* in Europe, it would be relevant to integrate future ABC
analysis of genotyped populations from the Italian peninsula, one of the three European glacial
refuges, and more oriental population.

Several recent studies have highlighted cryptic diversity within the Obsoletus/Scoticus
Complex (*C. obsoletus* clade O2, *C. obsoletus* clade dark, *C. montanus*) and raised questions
about the taxonomic status of *C. montanus*, a species phylogenetically close to *C. obsoletus*.
The taxonomic status of these clades is not yet fully established although there are strong
indications that these groups constitute species in their own right (Ander *et al*, 2013;
Meiswinkel *et al*, 2007; Mignotte *et al*, 2020). However, most of these studies were conducted
using mitochondrial markers. These markers, by responding differently than nuclear markers
to demographic fluctuations may have diverged more strongly while nuclear genes may have
retained variability. It is thus possible that several highly divergent mtDNA lineages coexist
over large geographical range despite the absence of complete reproductive isolation (Despres,
2019). Therefore, it is possible that by concentrating our study on *C. obsoletus sensu stricto*,
we may have missed part of the evolutionary history of this species.

The haplotype mismatch distribution, the lack of matrilineal subdivision, the star-like pattern
of median-joining haplotype network centred around few frequent haplotypes combined to
significant negative values of Fu's and Tajima's D suggest a rapid demographic expansion of
*C. obsoletus* in Europe. Recent demographic expansion is consistent with demographic history
of most European species modified during the glacial episodes of the Pleistocene (Bowler and
Benton, 2005; Hewitt, 2004; Phillips *et al*, 2010; Widmer and Lexer, 2001) but inconsistent
with our ABC results which report more recent events. A demographic expansion has been also
been detected in *C. imicola* with an increase of the effective population size in the western
Mediterranean basin and sub-Saharan Africa during the last 17,000 years and 80,000 years
respectively (Jacquet *et al*, 2015). Given the dependence of *Culicoides* on bovine hosts, the
development of modern human agriculture may have played an important role in the population
explosion of *C. obsoletus* in Europe. However, the domestication of ungulates has only taken
place in the Middle East and Europe during the last 10,000 years. It thus appears important to
be able to specify the dating of these expansions in order to be able to refine our hypotheses.

**Conclusion**

Our study highlights important gene flows at the continental level and the connectivity between
populations separated by a high geographical distance, possible thanks to the stratified
dispersion of *Culicoides* combining short and long-distance displacements. The use of
mitochondrial and microsatellite markers has allowed us to highlight a lack of phylogeographic
and genetic structure, a recent increase in effective population size, a recent divergence of
Western European *Culicoides* populations (5,000 to 10,000 years ago), and post-glacial origin
located in the Balkans and the arctic region of ancestral populations. The population expansion
highlights the potential role of the glacial episode and the post-glacial colonization event in
shaping the geographic distribution of *C. obsoletus*. The observed difference in structure
between nuclear and mitochondrial markers suggests that a coupled effect of stratified
*Culicoides* dispersal, displacement by human activities, and the recent increase in reproductive
activity may have led to genetic homogenization, erasing the past phylogeographic structure of
the population. However, a more precise time of expansion and divergence would help to better
understand the genetic homogeneity of *C. obsoletus* in the eastern Palearctic region. Our result
on the demographic and evolutionary history of *C. obsoletus* leads us to reassess the importance
of this vector in the dispersal of the BTV.
## Author contributions

A.M., K.H. and C.G. designed the study. I.G., M.D. contributed the identification of *Culicoides*. A.M., L.G., have performed molecular biology manipulations. A.M. analysed the data. K.H. and C.G contributed to the manuscript firstly written by A.M. All authors read and commented the final manuscript version.

## Acknowledgement

This study was partially funded by the EU grant H2020-727393 PALE-Blu and by the VectorNet project (OC/EFSA/AHAW/2013/02-FWC1) funded by the European Centre for Disease Prevention and Control (ECDC) and the European Food Safety Authority (EFSA). The contents of this publication are the sole responsibility of the authors and do not necessarily reflect the views of the European Commission, of the ECDC or of the EFSA. The authors are grateful to all European partners who gave assistance in *Culicoides* sampling and shipping.
Figures and tables

<figure><img src="image_34.png" /><figcaption>**Figure 1.** Graphical representation of the 4 scenarios simulated using an approximate Bayesian computation (ABC) approach on microsatellite data. Details of all scenarios and parameters are shown in Supplementary Table S2. Red circles correspond to a bottleneck event. T1, T2 and T3 correspond to successive dates of divergence times. NO: population from Norway (Narvik), SP: population from Spain (Torre del campo), FR: population from France (Les Hitres), GR: population from Greece (Vestino).</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>**Figure 2.** Median-joining haplotype networks of cox1 mitochondrial sequences of *C. obsoletus*. The size of the circles is proportional to the number of individuals with that haplotype. Colors represent the country of sampled individuals. The length of the branches separating haplotypes is proportional to the number of mutational steps between them.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>**Figure 3.** Mitochondrial genetic differentiation and isolation-by-distance analyses of *C. obsoletus*. Results of the smoothing of pairwise measures performed with MAPI and based on (a) FST/(1-FST) and (b) Isolation-by-distance analysis, Mantel test and a kernel density estimate of pairwise correlations from each sample between (FST/(1-FST) and geographic distance In panels (a), genetic dissimilarity is represented by a color scale ranging from brown (lower genetic dissimilarity) to blue (higher genetic dissimilarity). The Canary Islands populations are not considered in these analyses.</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>Figure 4. Mismatch distribution plots for 106 haplotypes of cox1 of *C. obsoletus*. (a) The mean estimate and the 90, 95 and 99% highest posterior density (HPD) limits are indicated.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>**Figure 5.** Microsatellite genetic clustering and genetic diversity of *C. obsoletus* sampled populations. Results of the genetic clustering analyses performed with STRUCTURE (a) and rarefied allelic richness of populations (b). Panels (1) corresponds to an enlargement of the geographical area (a) of the panel circled in red. One specific color has been assigned to each genetic cluster in panels (a) and (1).</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>**Figure 6.** Microsatellite genetic differentiation and isolation-by-distance analyses of *C. obsoletus*. Results of the smoothing of pairwise measures performed with MAPI and based on (a) FST/(1-FST) and (b) Isolation-by-distance analysis, Mantel test and a kernel density estimate of pairwise correlations from each sample between (FST/(1-FST) and geographic distance. In panels (a), genetic dissimilarity is represented by a color scale ranging from brown (lower genetic dissimilarity) to blue (higher genetic dissimilarity). The Canary Islands population is not considered in these analyses because of its insularity.</figcaption></figure>
<table><caption>Tableau 1. Result of confidence in scenario selection based on posterior probabilities to retrace evolutionary history of *C. obsoletus* in Europe by approximate Bayesian Computation (ABC) on microsatellite data. 95% confidence intervals and type I and II errors were presented. Red circles correspond to a bottleneck event. T1, T2 and T3 correspond to successive dates of divergence between two taxa. NO: Norwegian population (Narvik), SP: Spanish population (Torre del campo), FR: French population (Les Hitres), GR: Greek population (Velestino).</caption><thead><tr><th>Scenario</th><th>Description of tested scenarios</th><th>Posterior probability [95% CI]</th><th>Error of scenario choiceType I</th><th>Error of scenario choice Type II</th></tr></thead><tbody><tr><td>1</td><td></td><td>0.6350 [0.5795,0.6905]</td><td>0.432</td><td>0.142</td></tr><tr><td>2</td><td></td><td>0.2672 [0.2146,0.3198]</td><td></td><td></td></tr><tr><td>3</td><td></td><td>0.0033 [0.0000,0.0230]</td><td></td><td></td></tr><tr><td>4</td><td></td><td>0.0002 [0.0000,0.0201]</td><td></td><td></td></tr></tbody></table>
**Supplementary information**

**Tableau S1.** Metadata of the 2409 sequences of cox1 *C. obsoletus*.

The table being too voluminous it is not included in this manuscript.
<table><caption>**Tableau S2. Result of most probable scenario for all ABC simulation.** 95% confidence intervals and type I and II errors were presented to. Red circles corresponds to a bottleneck event. T1, T2 and T3 correspond to successive dates of divergence between two taxa. NO: Norway population (Narvik), SP: Spain population (Torre del campo), FR: French population (Les Hitres), GR: Greece populations (Vestino).</caption><thead><tr><th>Scenario</th><th>Description of selected scenarios</th><th>Divergence date</th><th>Posterior probability [95% CI]</th><th>Error of scenario choiceType I</th><th>Error of scenario choice Type II</th></tr></thead><tbody><tr><td>1</td><td></td><td>T1: 10 000-5000 years<br>T2: 5000-2500 years</td><td><strong>0.6350</strong><br><strong>[0.5795,0.6905]</strong></td><td>0.432</td><td>0.142</td></tr><tr><td>2</td><td></td><td>T1: 10 000-5000 years<br>T2: 5000-2500 years</td><td>0.2672<br><strong>[0.2146,0.3198]</strong></td><td></td><td></td></tr><tr><td>3</td><td></td><td>T1: 10 000-5000 years<br>T2: 5000-2500 years<br>T3: 2500-1 years</td><td>0.0033<br><strong>[0.0000,0.0230]</strong></td><td></td><td></td></tr><tr><td>4</td><td></td><td>T1: 10 000-5000 years<br>T2: 5000-2500 years<br>T3: 2500-1 years</td><td>0.0002<br><strong>[0.0000,0.0201]</strong></td><td></td><td></td></tr><tr><td>5</td><td></td><td>T1: 25 000-12 500 years<br>T2: 12 500-7500 years</td><td>0.0587<br><strong>[0.0391,0.0783]</strong></td><td></td><td></td></tr><tr><td>6</td><td></td><td>T1: 25 000-12 500 years<br>T2: 12 500-7500 years</td><td>0.0348<br><strong>[0.0170,0.0527]</strong></td><td></td><td></td></tr><tr><td>7</td><td></td><td>T1: 25 000-12 500 years<br>T2: 12 500-7500 years<br>T3:7500-5000 years</td><td>0.0004<br><strong>[0.0000,0.0203]</strong></td><td></td><td></td></tr><tr><td>8</td><td></td><td>T1: 25 000-12 500 years<br>T2: 12 500-7500 years<br>T3:7500-5000 years</td><td>0.0000<br><strong>[0.0000,0.0200]</strong></td><td></td><td></td></tr><tr><td>9</td><td></td><td>T1: 10 000-5000 years<br>T2: 5000-2500 years</td><td>0.0000<br><strong>[0.0000,0.0200]</strong></td><td></td><td></td></tr><tr><td>10</td><td></td><td>T1: 100 000-50 000 years<br>T2: 50 000-25 000 years</td><td>0.0000<br><strong>[0.0000,0.0200]</strong></td><td></td><td></td></tr><tr><td>11</td><td></td><td>T1: 100 000-50 000 years<br>T2: 50 000-25 000 years<br>T3: 25 000-20 000 years</td><td>0.0000<br><strong>[0.0000,0.0200]</strong></td><td></td><td></td></tr><tr><td>12</td><td></td><td>T1: 100 000-50 000 years<br>T2: 50 000-25 000 years<br>T3: 25 000-20 000 years</td><td>0.0000<br><strong>[0.0000,0.0200]</strong></td><td></td><td></td></tr></tbody></table>
Tableau S3. Demographic, historical and mutational prior parameters for ABC analyses.

Effective population sizes (N) are expressed in number of diploid individuals and times of events (t) in number of generations going back to the past.

<table><thead><tr><th>Parameter name</th><th>Distribution Prior</th><th>Parameters abbreviation</th></tr></thead><tbody><tr><td>Mean mutation rate for nucleotide instability</td><td>Uniform [10-8 ; 10-4 ]</td><td>SNI</td></tr><tr><td>Mean parameter of the geometric distribution of length</td><td>Uniform [10-1 ; 3.10-1 ]</td><td>P</td></tr><tr><td>Mean microsatellite mutation model</td><td>Loguniform [6.10-6 ; 10-4 ]</td><td>μ</td></tr><tr><td>Very recent ancestral divergence time</td><td>Uniform [40 000; 20 000]</td><td>T1</td></tr><tr><td>Time of second event of divergence</td><td>Uniform [20 000; 10 000]</td><td>T2</td></tr><tr><td>Time of third event of divergence</td><td>Uniform [10 000; 4]</td><td>T3</td></tr><tr><td>Recent ancestral divergence time</td><td>Uniform [100 000; 50 000]</td><td>T1</td></tr><tr><td>Time of second event of divergence</td><td>Uniform [50 000; 30 000]</td><td>T2</td></tr><tr><td>Time of third event of divergence</td><td>Uniform [30 000; 20 000]</td><td>T3</td></tr><tr><td>Old ancestral divergence time</td><td>Uniform [400 000; 200 000]</td><td>T1</td></tr><tr><td>Time of second event of divergence</td><td>Uniform [200 000; 100 000]</td><td>T2</td></tr><tr><td>Time of third event of divergence</td><td>Uniform [100 000; 80 000]</td><td>T3</td></tr><tr><td>Duration of bottleneck</td><td>Uniform [0; 100]</td><td>DB</td></tr><tr><td>Effective population size during bottleneck</td><td>Loguniform [2; 1000]</td><td>Nbott</td></tr><tr><td>Stable effective population size in ancestral (a) and European (e) populations</td><td>Loguniform [100; 200000]</td><td>Na, Ne</td></tr></tbody></table>
Tableau S4. Population genetic diversity of 2409 cox1 sequence of *C. obsoletus* in Europe.

Pi: Nucleotide diversity per site, Nh: haplotype number.

<table><thead><tr><th>Locality</th><th>Country</th><th>Longitude</th><th>Latitude</th><th>Collection date</th><th>Abbreviation</th><th>Nseq</th><th>Number of variable sites</th><th>Nh</th><th>Haplotype diversity</th><th>Standard Deviation of Haplotype diversity</th><th>Pi</th><th>Standard deviation of Pi</th><th>Ar</th><th>Fu's Fs statistic</th><th>Tajima's D</th></tr></thead><tbody><tr><td>Montana</td><td>Bulgaria</td><td>23,23</td><td>43,41</td><td>2016</td><td>Mon</td><td>143</td><td>22</td><td>21</td><td>0,737</td><td>0,034</td><td>0,004</td><td>0,000</td><td>2,846</td><td>-7,127</td><td>-1,100</td></tr><tr><td>Nikolovo</td><td>Bulgaria</td><td>26,115334</td><td>43,87123</td><td>2016</td><td>Nik</td><td>146</td><td>20</td><td>18</td><td>0,697</td><td>0,037</td><td>0,004</td><td>0,000</td><td>2,849</td><td>-5,454</td><td>-1,147</td></tr><tr><td>Dyssevej</td><td>Denmark</td><td>12,11</td><td>55,2</td><td>2015</td><td>Dys</td><td>83</td><td>30</td><td>24</td><td>0,873</td><td>0,026</td><td>0,007</td><td>0,000</td><td>2,036</td><td>-7,930</td><td>-1,113</td></tr><tr><td>Grøddevej</td><td>Denmark</td><td>9,18</td><td>56,13</td><td>2015</td><td>Grw</td><td>50</td><td>22</td><td>17</td><td>0,759</td><td>0,062</td><td>0,005</td><td>0,001</td><td>2,540</td><td>-5,983</td><td>-1,308</td></tr><tr><td>Leestrupvej</td><td>Denmark</td><td>12,05</td><td>55,22</td><td>2015</td><td>Lee</td><td>92</td><td>23</td><td>20</td><td>0,752</td><td>0,042</td><td>0,005</td><td>0,001</td><td>2,348</td><td>-5,543</td><td>-1,124</td></tr><tr><td>Lægdsgårdsvej</td><td>Denmark</td><td>9,12</td><td>56,12</td><td>2015</td><td>Leg</td><td>108</td><td>30</td><td>26</td><td>0,823</td><td>0,029</td><td>0,006</td><td>0,000</td><td>2,259</td><td>-10,434</td><td>-1,338</td></tr><tr><td>Inari</td><td>Finland</td><td>27,56</td><td>68,67</td><td>2015</td><td>Ina</td><td>8</td><td>3,000</td><td>3</td><td>0,750</td><td>0,096</td><td>0,003</td><td>0,000</td><td>5,750</td><td>0,971</td><td>1,220</td></tr><tr><td>Sotkamo</td><td>Finland</td><td>28,72</td><td>64,04</td><td>2015</td><td>Sot</td><td>19</td><td>9</td><td>5</td><td>0,737</td><td>0,069</td><td>0,007</td><td>0,001</td><td>3,211</td><td>2,480</td><td>0,962</td></tr><tr><td>Barjols</td><td>France</td><td>5,98</td><td>43,54</td><td>2016</td><td>Bar</td><td>13</td><td>5</td><td>5</td><td>0,859</td><td>0,050</td><td>0,003</td><td>0,000</td><td>2,000</td><td>0,043</td><td>0,298</td></tr><tr><td>Bonnay</td><td>France</td><td>2,51</td><td>49,93</td><td>2016</td><td>Bon</td><td>38</td><td>16</td><td>12</td><td>0,774</td><td>0,056</td><td>0,005</td><td>0,001</td><td>2,368</td><td>-2,711</td><td>-0,992</td></tr><tr><td>Les Hitres</td><td>France</td><td>4,07</td><td>45,75</td><td>2016</td><td>Les</td><td>43</td><td>11</td><td>10</td><td>0,759</td><td>0,053</td><td>0,003</td><td>0,000</td><td>2,698</td><td>-2,223</td><td>-0,830</td></tr><tr><td>Reygnac</td><td>France</td><td>1,8176424</td><td>45,3429047</td><td>2016</td><td>Rey</td><td>11</td><td>14</td><td>6</td><td>0,800</td><td>0,114</td><td>0,008</td><td>0,002</td><td>1,909</td><td>0,494</td><td>-0,483</td></tr><tr><td>Sucé-sur-Erdre</td><td>France</td><td>-1,54</td><td>47,32</td><td>2016</td><td>Suc</td><td>8</td><td>8</td><td>4</td><td>0,750</td><td>0,139</td><td>0,005</td><td>0,002</td><td>2,125</td><td>0,896</td><td>-0,866</td></tr><tr><td>Peetzig</td><td>Germany</td><td>13,94</td><td>53,08</td><td>2016</td><td>Pee</td><td>8</td><td>3</td><td>3</td><td>0,607</td><td>0,164</td><td>0,002</td><td>0,001</td><td>3,375</td><td>0,723</td><td>0,585</td></tr><tr><td>Schleiz</td><td>Germany</td><td>11,81</td><td>50,57</td><td>2016</td><td>Sch</td><td>33</td><td>11</td><td>8</td><td>0,680</td><td>0,054</td><td>0,003</td><td>0,001</td><td>3,121</td><td>-1,655</td><td>-1,349</td></tr><tr><td>Ziegelbach</td><td>Germany</td><td>9,86</td><td>47,89</td><td>2016</td><td>Zie</td><td>31</td><td>9</td><td>5</td><td>0,652</td><td>0,063</td><td>0,004</td><td>0,001</td><td>3,194</td><td>2,018</td><td>-0,116</td></tr><tr><td>Abercorn</td><td>Great Britain</td><td>-3,455284</td><td>55,995067</td><td>2014</td><td>Abe</td><td>24</td><td>11</td><td>7</td><td>0,779</td><td>0,057</td><td>0,004</td><td>0,001</td><td>2,625</td><td>0,073</td><td>-0,553</td></tr><tr><td>Blackmoor</td><td>Great Britain</td><td>-3,94</td><td>51,18</td><td>2015</td><td>Bla</td><td>30</td><td>13</td><td>7</td><td>0,756</td><td>0,050</td><td>0,005</td><td>0,001</td><td>2,567</td><td>0,992</td><td>-0,400</td></tr><tr><td>London</td><td>Great Britain</td><td>-0,13</td><td>51,51</td><td>2010</td><td>Lon</td><td>27</td><td>16</td><td>8</td><td>0,604</td><td>0,107</td><td>0,005</td><td>0,001</td><td>4,778</td><td>-0,452</td><td>-1,328</td></tr><tr><td>Thorpe</td><td>Great Britain</td><td>-1,78</td><td>53,04</td><td>2015</td><td>Tho</td><td>17</td><td>15</td><td>9</td><td>0,831</td><td>0,085</td><td>0,006</td><td>0,001</td><td>2,471</td><td>-1,744</td><td>-0,992</td></tr><tr><td>Varybobi</td><td>Greece</td><td>23,78</td><td>38,13</td><td>2016</td><td>Var</td><td>50</td><td>20</td><td>19</td><td>0,863</td><td>0,035</td><td>0,006</td><td>0,001</td><td>3,000</td><td>-9,043</td><td>-1,242</td></tr><tr><td>Velestino</td><td>Greece</td><td>22,75</td><td>39,38</td><td>2016</td><td>Vel</td><td>112</td><td>19</td><td>17</td><td>0,731</td><td>0,035</td><td>0,004</td><td>0,000</td><td>4,536</td><td>-5,840</td><td>-1,448</td></tr><tr><td>Brunico</td><td>Italy</td><td>11,96</td><td>46,8</td><td>2016</td><td>Bru</td><td>46</td><td>20</td><td>12</td><td>0,850</td><td>0,029</td><td>0,006</td><td>0,001</td><td>1,957</td><td>-1,095</td><td>-0,813</td></tr><tr><td>Mammola</td><td>Italy</td><td>16,19</td><td>38,38</td><td>2016</td><td>Mam</td><td>19</td><td>10</td><td>6</td><td>0,760</td><td>0,070</td><td>0,007</td><td>0,001</td><td>1,842</td><td>1,917</td><td>1,255</td></tr><tr><td>Teramo</td><td>Italy</td><td>13,67</td><td>42,64</td><td>2016</td><td>Ter</td><td>14</td><td>5</td><td>4</td><td>0,571</td><td>0,132</td><td>0,003</td><td>0,001</td><td>3,143</td><td>0,531</td><td>-0,099</td></tr><tr><td>Kalnciemi</td><td>Latvia</td><td>25,44</td><td>57,66</td><td>2016</td><td>Kal</td><td>165</td><td>13</td><td>11</td><td>0,431</td><td>0,044</td><td>0,003</td><td>0,000</td><td>4,691</td><td>-0,284</td><td>-0,840</td></tr><tr><td>Mazlimbas</td><td>Latvia</td><td>25,4</td><td>57,71</td><td>2016</td><td>Maz</td><td>75</td><td>11</td><td>7</td><td>0,646</td><td>0,042</td><td>0,006</td><td>0,000</td><td>4,627</td><td>2,692</td><td>0,645</td></tr><tr><td>Oubogo</td><td>Macedonia</td><td>21,99037</td><td>41,625636</td><td>2016</td><td>Oub</td><td>34</td><td>18</td><td>12</td><td>0,832</td><td>0,054</td><td>0,005</td><td>0,001</td><td>2,206</td><td>-2,931</td><td>-1,272</td></tr><tr><td>Bab Taza</td><td>Morocco</td><td>-5,189793</td><td>31,264597</td><td>2016</td><td>Bab</td><td>1</td><td></td><td>1</td><td></td><td></td><td></td><td></td><td>5,000</td><td></td><td></td></tr><tr><td>Chefchaouen</td><td>Morocco</td><td>-5,273638</td><td>35,161626</td><td>2016</td><td>Che</td><td>16</td><td>8</td><td>4</td><td>0,625</td><td>0,093</td><td>0,006</td><td>0,001</td><td>2,688</td><td>4,400</td><td>0,993</td></tr><tr><td>Dardara</td><td>Morocco</td><td>-5,301693</td><td>34,489081</td><td>2016</td><td>Dar</td><td>2</td><td>2</td><td>2</td><td>1,000</td><td>0,500</td><td>0,003</td><td>0,002</td><td>2,000</td><td>0,693</td><td></td></tr><tr><td>Moulay Bousselham</td><td>Morocco</td><td>-6,245313</td><td>31,264597</td><td>2016</td><td>Mou</td><td>1</td><td></td><td>1</td><td></td><td></td><td></td><td></td><td>1,000</td><td></td><td></td></tr><tr><td>Oulad Mtaa</td><td>Morocco</td><td>-8,256511</td><td>31,264597</td><td>2016</td><td>Oul</td><td>1</td><td></td><td>1</td><td></td><td></td><td></td><td></td><td>1,000</td><td></td><td></td></tr><tr><td>Rabat</td><td>Morocco</td><td>-6,8</td><td>33,96</td><td>2010</td><td>Rab</td><td>4</td><td>4</td><td>3</td><td>0,833</td><td>0,222</td><td>0,004</td><td>0,001</td><td>1,000</td><td>0,461</td><td>1,365</td></tr><tr><td>Tizgzaouine</td><td>Morocco</td><td>-8,130402</td><td>31,767321</td><td>2016</td><td>Tiz</td><td>1</td><td></td><td>1</td><td></td><td></td><td></td><td></td><td>3,000</td><td></td><td></td></tr><tr><td>Lomm</td><td>Netherlands</td><td>6,2</td><td>51,45</td><td>2017</td><td>Lom</td><td>18</td><td>12</td><td>8</td><td>0,797</td><td>0,088</td><td>0,004</td><td>0,001</td><td>2,444</td><td>-1,860</td><td>-1,363</td></tr><tr><td>Olst</td><td>Netherlands</td><td>6,12</td><td>52,37</td><td>2017</td><td>Ols</td><td>11</td><td>11</td><td>5</td><td>0,709</td><td>0,137</td><td>0,004</td><td>0,002</td><td>2,909</td><td>0,315</td><td>-1,471</td></tr><tr><td>Narvik</td><td>Norway</td><td>17,51</td><td>68,48</td><td>2015</td><td>Nar</td><td>82</td><td>13</td><td>8</td><td>0,684</td><td>0,029</td><td>0,003</td><td>0,000</td><td>5,098</td><td>-0,180</td><td>-1,042</td></tr><tr><td>Storforshei</td><td>Norway</td><td>14,81</td><td>66,31</td><td>2015</td><td>Sto</td><td>4</td><td>0</td><td>1</td><td></td><td></td><td></td><td></td><td>2,750</td><td></td><td></td></tr><tr><td>Tiller</td><td>Norway</td><td>10,42</td><td>63,34</td><td>2015</td><td>Til</td><td>19</td><td>3</td><td>4</td><td>0,626</td><td>0,067</td><td>0,002</td><td>0,000</td><td>6,053</td><td>0,462</td><td>0,275</td></tr><tr><td>Jaksamanice</td><td>Poland</td><td>22,86</td><td>49,751</td><td>2016</td><td>Jak</td><td>93</td><td>15</td><td>13</td><td>0,678</td><td>0,035</td><td>0,003</td><td>0,000</td><td>2,882</td><td>-4,184</td><td>-1,326</td></tr><tr><td>Olszewnica</td><td>Poland</td><td>22,51</td><td>51,76</td><td>2016</td><td>Olz</td><td>25</td><td>8</td><td>6</td><td>0,690</td><td>0,065</td><td>0,003</td><td>0,001</td><td>3,400</td><td>0,671</td><td>-0,849</td></tr><tr><td>Pikulice</td><td>Poland</td><td>22,79</td><td>49,75</td><td>2016</td><td>Pik</td><td>195</td><td>18</td><td>15</td><td>0,671</td><td>0,024</td><td>0,003</td><td>0,000</td><td>3,815</td><td>-4,611</td><td>-1,347</td></tr></tbody></table>
<table><tr><td>Aboadela</td><td>Portugal</td><td>-8</td><td>41,27</td><td>2009</td><td>Abo</td><td>26</td><td>14</td><td>12</td><td>0,880</td><td>0,034</td><td>0,005</td><td>0,001</td><td>1,962</td><td>-2,211</td><td>-1,048</td></tr><tr><td>Bajouca</td><td>Portugal</td><td>-8,77</td><td>39,87</td><td>2009</td><td>Baj</td><td>1</td><td></td><td>1</td><td></td><td></td><td></td><td></td><td>1,000</td><td></td><td></td></tr><tr><td>Aleksandrovo</td><td>Serbia</td><td>20,58</td><td>45,63</td><td>2016</td><td>Ale</td><td>61</td><td>12</td><td>9</td><td>0,649</td><td>0,041</td><td>0,004</td><td>0,000</td><td>4,525</td><td>-0,884</td><td>-0,834</td></tr><tr><td>Cirikovac</td><td>Serbia</td><td>21,19</td><td>44,67</td><td>2016</td><td>Cir</td><td>78</td><td>12</td><td>12</td><td>0,680</td><td>0,034</td><td>0,003</td><td>0,000</td><td>2,423</td><td>-2,552</td><td>-0,639</td></tr><tr><td>Gamzigrad</td><td>Serbia</td><td>22,18</td><td>43,91</td><td>2016</td><td>Gam</td><td>71</td><td>14</td><td>12</td><td>0,743</td><td>0,036</td><td>0,004</td><td>0,000</td><td>2,732</td><td>-2,199</td><td>-0,780</td></tr><tr><td>Pacir</td><td>Serbia</td><td>19,45</td><td>45,89</td><td>2016</td><td>Pac</td><td>59</td><td>12</td><td>6</td><td>0,704</td><td>0,036</td><td>0,005</td><td>0,001</td><td>5,729</td><td>2,542</td><td>-0,153</td></tr><tr><td>Arauut</td><td>Spain</td><td>1,760802</td><td>41,965161</td><td>2017</td><td>Ara</td><td>30</td><td>12</td><td>14</td><td>0,864</td><td>0,048</td><td>0,004</td><td>0,001</td><td>1,700</td><td>-4,367</td><td>-0,869</td></tr><tr><td>Arucas</td><td>Spain</td><td>-15,5</td><td>28,13</td><td>2010</td><td>Aru</td><td>37</td><td>2</td><td>3</td><td>0,107</td><td>0,068</td><td>0,000</td><td>0,000</td><td>4,054</td><td>-2,626</td><td>-1,493</td></tr><tr><td>Castril</td><td>Spain</td><td>-2,34</td><td>37,73</td><td>2017</td><td>Cas</td><td>16</td><td>9</td><td>5</td><td>0,683</td><td>0,091</td><td>0,003</td><td>0,001</td><td>3,063</td><td>0,237</td><td>-1,327</td></tr><tr><td>Orgiva</td><td>Spain</td><td>-2,99</td><td>36,85</td><td>2017</td><td>Org</td><td>6</td><td>2</td><td>2</td><td>0,600</td><td>0,129</td><td>0,002</td><td>0,000</td><td>3,667</td><td>1,938</td><td>1,753</td></tr><tr><td>Palma Torre del Campo</td><td>Spain</td><td>2,65</td><td>39,56</td><td>2015</td><td>Pal</td><td>26</td><td>9</td><td>8</td><td>0,723</td><td>0,062</td><td>0,004</td><td>0,001</td><td>3,308</td><td>-1,469</td><td>-0,694</td></tr><tr><td>Dittingen</td><td>Switzerland</td><td>7,5</td><td>47,44</td><td>2009</td><td>Dit</td><td>13</td><td>8</td><td>4</td><td>0,731</td><td>0,079</td><td>0,005</td><td>0,001</td><td>2,923</td><td>1,948</td><td>-0,161</td></tr><tr><td>Muhleturnen</td><td>Switzerland</td><td>7,51</td><td>46,81</td><td>2009</td><td>Muh</td><td>18</td><td>13</td><td>9</td><td>0,895</td><td>0,048</td><td>0,005</td><td>0,001</td><td>2,056</td><td>-1,904</td><td>-0,832</td></tr><tr><td>Şahinköy</td><td>Turkey</td><td>28,237</td><td>40,235</td><td>2016</td><td>gah</td><td>15</td><td>12</td><td>5</td><td>0,676</td><td>0,105</td><td>0,004</td><td>0,001</td><td>3,333</td><td>0,894</td><td>-1,366</td></tr><tr><td>Kuzgöl</td><td>Turkey</td><td>32,09938</td><td>40,91156</td><td>2016</td><td>kuz</td><td>31</td><td>23</td><td>12</td><td>0,882</td><td>0,039</td><td>0,007</td><td>0,001</td><td>1,355</td><td>-1,650</td><td>-1,116</td></tr><tr><td>All data Set</td><td></td><td></td><td></td><td></td><td></td><td>2409</td><td>67</td><td>106</td><td>0,772</td><td>0,005</td><td>0,005</td><td></td><td>1,000</td><td>-30,374</td><td>-1,985</td></tr></table>
<table><caption>Tableau S5. Primers of the 11 microsatellite markers used to genotype *C. obsoletus* populations. “DYE”: fluorochrome name, “Tm”: half denaturation temperature in degree, “bp”: base pairs.</caption><thead><tr><th>Multiplex</th><th>Locus</th><th>DYE</th><th>Repeat motif</th><th>Forward primer (5'-3')</th><th>Tm</th><th>Reverse primer (3'-5')</th><th>Tm</th><th>Size range (bp)</th><th>Accession numbers</th></tr></thead><tbody><tr><td rowspan="5">Plex1</td><td>OBSms5</td><td>FAM</td><td>(TG)15</td><td>GCTCATTTATGTAATTTGTGCC</td><td>56,4</td><td>CGAGACGGCATGTTTATCT</td><td>56,3</td><td>140-204</td><td>MT828835</td></tr><tr><td>OBSms11</td><td>FAM</td><td>(CT)12</td><td>CATCCATCCACTTTGCTCT</td><td>56,1</td><td>GGACTGTCATTCAAACAAACA</td><td>56,1</td><td>262-293</td><td>MT828838</td></tr><tr><td>OBSms2</td><td>NED</td><td>(TG)14</td><td>TGATAACAGATAACACGGAAGG</td><td>56,8</td><td>TGATGCAATAATGTCCTCGT</td><td>58,6</td><td>135-171</td><td>MT828832</td></tr><tr><td>OBSms28</td><td>NED</td><td>(CA)8</td><td>GCCGCAATTTGTGTTTT</td><td>55,5</td><td>TTCCTTGACCTGCGAAC</td><td>55,5</td><td>239-264</td><td>MT828842</td></tr><tr><td>OBSms3</td><td>VIC</td><td>(CA)13</td><td>ATTTACGCGACGATCATTT</td><td>55,8</td><td>AAAGCTCGAATACAGCACAA</td><td>56,2</td><td>140-184</td><td>MT828833</td></tr><tr><td>OBSms4</td><td>PET</td><td>(GA)12</td><td>AATGTAGGCCAAACAGAAA</td><td>53,3</td><td>TTGCATAACAACATCATCAC</td><td>52,7</td><td>160-180</td><td>MT828834</td></tr><tr><td rowspan="4">Plex2</td><td>OBSms29</td><td>FAM</td><td>(CA)5</td><td>GTTGACCGGGACACTTTT</td><td>55,8</td><td>CATGACAAGCACGAACG</td><td>55,4</td><td>252-268</td><td>MT828843</td></tr><tr><td>OBSms6</td><td>NED</td><td>(TA)12</td><td>GACCTAATTCATAATCGTCCAA</td><td>55,8</td><td>TTAAAGCACGTTCGATTTGT</td><td>56,1</td><td>175-201</td><td>MT828836</td></tr><tr><td>OBSms13</td><td>VIC</td><td>(CA)12</td><td>CGTTCCATCGACTCATCTT</td><td>56,1</td><td>GCTTTTACTCCTTCGTGATGT</td><td>56,1</td><td>249-287</td><td>MT828839</td></tr><tr><td>OBSms14</td><td>VIC</td><td>(AT)8</td><td>TCAAATCATTCCAAGAGCA</td><td>55,1</td><td>CCTCTCTGTGGTCTTTCGT</td><td>55,2</td><td>139-164</td><td>MT828840</td></tr><tr><td>OBSms10</td><td>PET</td><td>(AC)17</td><td>GACATGGAATCCTTAATTGG</td><td>54,5</td><td>GATTGCTTCAGACAAACACA</td><td>54,9</td><td>224-260</td><td>MT828837</td></tr></tbody></table>
**Tableau S6.** Microsatellite dataset of 829 *C. obsoletus* in Europe.

The table being too voluminous it is not included in this manuscript.
<table><caption>Tableau S7. Sampling sites and associated genetic diversity for microsatellite data set. N: sample size, Ho: observed heterozygosity, Hs: expected heterozygosity, FIS: Fixation index, Ar: allelic richness, Ne: sample size. P-values obtained for the Hardy-Weinberg test were adjusted for multiple comparisons at the nominal level (5%) of 0.00017.</caption><thead><tr><th>Population name</th><th>Country</th><th>Longitude</th><th>Latitude</th><th>N</th><th>Ho</th><th>Hs</th><th>Fis</th><th>Rarefied Ar</th><th>Ne (0.05)</th><th>95% CIs for Ne (Parametric)</th><th>95% CIs for Ne (JackKnife on Samples)</th></tr></thead><tbody><tr><td>Les Hitres</td><td>France</td><td>4,070</td><td>45,750</td><td>35</td><td>0,729</td><td>0,796</td><td>0,085</td><td>7,792</td><td>196.3</td><td>82,0 - Infinite</td><td>60,8 - Infinite</td></tr><tr><td>Bonnay</td><td>France</td><td>2,510</td><td>49,930</td><td>30</td><td>0,703</td><td>0,800</td><td>0,121</td><td>7,647</td><td>263.4</td><td>83,4 - Infinite</td><td>56,9 - Infinite</td></tr><tr><td>Nikolovo</td><td>Bulgaria</td><td>26,115</td><td>43,871</td><td>29</td><td>0,752</td><td>0,790</td><td>0,047</td><td>7,648</td><td>75.2</td><td>44,8 - 189,1</td><td>32,3 - Infinite</td></tr><tr><td>Montana</td><td>Bulgaria</td><td>23,230</td><td>43,410</td><td>20</td><td>0,750</td><td>0,777</td><td>0,035</td><td>6,873</td><td>Infinite</td><td>76,5 - Infinite</td><td>58,0 - Infinite</td></tr><tr><td>Arucas</td><td>Spain</td><td>-15,500</td><td>28,130</td><td>34</td><td>0,722</td><td>0,743</td><td>0,029</td><td>5,903</td><td>89.9</td><td>48,3 - 342,2</td><td>37,2 - Infinite</td></tr><tr><td>Dyssevej</td><td>Denmark</td><td>12,110</td><td>55,200</td><td>31</td><td>0,718</td><td>0,777</td><td>0,075</td><td>7,464</td><td>256.6</td><td>73,7 - Infinite</td><td>51,5 - Infinite</td></tr><tr><td>Grøddevej</td><td>Denmark</td><td>9,180</td><td>56,130</td><td>38</td><td>0,740</td><td>0,771</td><td>0,041</td><td>7,359</td><td>183.2</td><td>80,0 - Infinite</td><td>80,0 - Infinite</td></tr><tr><td>Blackmoor Torre del Campo</td><td>Great Britain</td><td>-3,940</td><td>51,180</td><td>23</td><td>0,688</td><td>0,773</td><td>0,109</td><td>7,385</td><td>25.2</td><td>16,9 - 43,1</td><td>9,6 - 1368,4</td></tr><tr><td>Ziegelbach</td><td>Germany</td><td>9,860</td><td>47,890</td><td>31</td><td>0,757</td><td>0,794</td><td>0,047</td><td>7,836</td><td>234.4</td><td>82,1 - Infinite</td><td>79,7 - Infinite</td></tr><tr><td>Velestino</td><td>Greece</td><td>22,750</td><td>39,380</td><td>30</td><td>0,718</td><td>0,819</td><td>0,124</td><td>7,546</td><td>42.7</td><td>29,6 - 69,7</td><td>24,4 - 107,9</td></tr><tr><td>Brunico</td><td>Italy</td><td>11,960</td><td>46,800</td><td>58</td><td>0,686</td><td>0,789</td><td>0,131</td><td>7,692</td><td>200.0</td><td>108,1 - 823,2</td><td>72,4 - Infinite</td></tr><tr><td>Mazlimbas</td><td>Latvia</td><td>25,400</td><td>57,710</td><td>40</td><td>0,747</td><td>0,797</td><td>0,062</td><td>7,811</td><td>80.7</td><td>52,7 - 152,3</td><td>34,6 - Infinite</td></tr><tr><td>Kalnciemi</td><td>Latvia</td><td>25,440</td><td>57,660</td><td>39</td><td>0,744</td><td>0,781</td><td>0,047</td><td>7,343</td><td>Infinite</td><td>362,9 - Infinite</td><td>211,9 - Infinite</td></tr><tr><td>Narvik</td><td>Norway</td><td>17,510</td><td>68,480</td><td>50</td><td>0,757</td><td>0,784</td><td>0,034</td><td>6,773</td><td>76.8</td><td>51,4 - 135,2 -</td><td>40,0 - 277,4</td></tr><tr><td>Aboadela</td><td>Portugal</td><td>-8,000</td><td>41,270</td><td>24</td><td>0,659</td><td>0,744</td><td>0,115</td><td>6,951</td><td>Infinite</td><td>85,6 - Infinite</td><td>77,8 - Infinite</td></tr><tr><td>Lomm</td><td>Netherlands</td><td>6,200</td><td>51,450</td><td>19</td><td>0,736</td><td>0,785</td><td>0,063</td><td>7,657</td><td>139.0</td><td>48,1 - Infinite</td><td>44,7 - Infinite</td></tr><tr><td>Pikulice</td><td>Poland</td><td>22,790</td><td>49,750</td><td>37</td><td>0,734</td><td>0,809</td><td>0,093</td><td>7,858</td><td>Infinite</td><td>250,1 - Infinite</td><td>187,1 - Infinite</td></tr><tr><td>Jaksamanice</td><td>Poland</td><td>22,860</td><td>49,751</td><td>32</td><td>0,730</td><td>0,799</td><td>0,086</td><td>7,662</td><td>Infinite</td><td>146,8 - Infinite</td><td>79,9 - Infinite</td></tr><tr><td>Cirikovac</td><td>Serbia</td><td>21,190</td><td>44,670</td><td>27</td><td>0,722</td><td>0,771</td><td>0,064</td><td>7,203</td><td>Infinite</td><td>98,9 - Infinite</td><td>78,9 - Infinite</td></tr><tr><td>Pacir</td><td>Serbia</td><td>19,450</td><td>45,890</td><td>32</td><td>0,728</td><td>0,791</td><td>0,08</td><td>7,671</td><td>82.9</td><td>45,6 - 276,9</td><td>33,1 - Infinite</td></tr><tr><td>Gamzigrad</td><td>Serbia</td><td>22,180</td><td>43,910</td><td>31</td><td>0,742</td><td>0,801</td><td>0,074</td><td>7,953</td><td>103.2</td><td>49,6 - 1263,9</td><td>44,5 - Infinite</td></tr><tr><td>Olszewnica</td><td>Poland</td><td>22,510</td><td>51,760</td><td>18</td><td>0,747</td><td>0,834</td><td>0,104</td><td>8,493</td><td>13492.6</td><td>75,9 - Infinite</td><td>48,4 - Infinite</td></tr><tr><td>Oubogo</td><td>Macedonia</td><td>21,990</td><td>41,626</td><td>28</td><td>0,747</td><td>0,824</td><td>0,093</td><td>8,449</td><td>73.7</td><td>43,0 - 199,1</td><td>33,8 - 2310,1</td></tr><tr><td>Abercorn</td><td>Great Britain</td><td>-3,466</td><td>55,997</td><td>21</td><td>0,723</td><td>0,809</td><td>0,107</td><td>7,590</td><td>918.7</td><td>53,9 - Infinite</td><td>36,4 - Infinite</td></tr><tr><td>Schleiz</td><td>Germany</td><td>11,810</td><td>50,570</td><td>31</td><td>0,736</td><td>0,807</td><td>0,088</td><td>7,649</td><td>190.3</td><td>69,7 - Infinite</td><td>61,5 - Infinite</td></tr><tr><td>Şahinköy</td><td>Turkey</td><td>28,237</td><td>40,235</td><td>25</td><td>0,716</td><td>0,829</td><td>0,136</td><td>8,513</td><td>123.5</td><td>55,7 - Infinite</td><td>39,2 - Infinite</td></tr></tbody></table>
<table><caption>**Tableau S8. Genetic diversity by locus.** Ho: observed heterozygosity, Hs: expected heterozygosity, FIS: Fixation index, SAD: short alleles dominance, CI: confidence index.</caption><thead><tr><th>Locus name</th><th>Alleles number</th><th>Ho</th><th>Hs</th><th>Fis</th><th>SAD r<sup>2</sup></th><th>SAD pvalue</th><th>Null alleles frequency [95% CI]</th></tr></thead><tbody><tr><td>OBSms11</td><td>29</td><td>0,810</td><td>0,872</td><td>0,071</td><td>0,375</td><td>0,054</td><td>0,043 [0,028-0,058]</td></tr><tr><td>OBSms2</td><td>29</td><td>0,781</td><td>0,849</td><td>0,080</td><td>-0,064</td><td>0,755</td><td>0,056 [0,041-0,074]</td></tr><tr><td>OBSms28</td><td>21</td><td>0,702</td><td>0,729</td><td>0,036</td><td>0,194</td><td>0,413</td><td>0,024 [0,004-0,046]</td></tr><tr><td>OBSms3</td><td>34</td><td>0,734</td><td>0,788</td><td>0,068</td><td>-0,322</td><td>0,068</td><td>0,051 [0,031-0,069]</td></tr><tr><td>OBSms4</td><td>20</td><td>0,778</td><td>0,816</td><td>0,047</td><td>0,111</td><td>0,661</td><td>0,031 [0,013-0,049]</td></tr><tr><td>OBSms5</td><td>47</td><td>0,783</td><td>0,890</td><td>0,120</td><td>-0,098</td><td>0,528</td><td>0,076 [0,058-0,094]</td></tr><tr><td>OBSms10</td><td>26</td><td>0,750</td><td>0,876</td><td>0,145</td><td>0,036</td><td>0,865</td><td>0,083 [0,063-0,102]</td></tr><tr><td>OBSms13</td><td>29</td><td>0,872</td><td>0,890</td><td>0,021</td><td>-0,256</td><td>0,198</td><td>0,018 [0,005-0,029]</td></tr><tr><td>OBSms14</td><td>19</td><td>0,700</td><td>0,760</td><td>0,079</td><td>-0,298</td><td>0,244</td><td>0,046 [0,025-0,066]</td></tr><tr><td>OBSms29</td><td>14</td><td>0,350</td><td>0,416</td><td>0,158</td><td>-0,081</td><td>0,794</td><td>0,114 [0,078-0,151]</td></tr><tr><td>OBSms6</td><td>40</td><td>0,724</td><td>0,827</td><td>0,125</td><td>-0,074</td><td>0,650</td><td>0,089 [0,069-0,109]</td></tr></tbody></table>
<table><caption>**Tableau S9.** Posterior checking using all the summary statistics for ABC inferences of scenario 1, the most probable scenarios.</caption><thead><tr><th>Summary statistics</th><th>Observed value</th><th>Proportion (simulated&lt;observed)</th></tr></thead><tbody><tr><td>NAL_1_1</td><td>10.3636</td><td>0.7180</td></tr><tr><td>NAL_1_2</td><td>12.0909</td><td>0.8390</td></tr><tr><td>NAL_1_3</td><td>11.2727</td><td>0.8200</td></tr><tr><td>NAL_1_4</td><td>8.2727</td><td>0.6945</td></tr><tr><td>HET_1_1</td><td>0.7835</td><td>0.6565</td></tr><tr><td>HET_1_2</td><td>0.7951</td><td>0.5030</td></tr><tr><td>HET_1_3</td><td>0.8178</td><td>0.7185</td></tr><tr><td>HET_1_4</td><td>0.7855</td><td>0.5625</td></tr><tr><td>VAR_1_1</td><td>10.9102</td><td>0.8370</td></tr><tr><td>VAR_1_2</td><td>12.1763</td><td>0.8530</td></tr><tr><td>VAR_1_3</td><td>8.6675</td><td>0.7030</td></tr><tr><td>VAR_1_4</td><td>9.4842</td><td>0.7630</td></tr><tr><td>MGW_1_1</td><td>0.8261</td><td>0.2085</td></tr><tr><td>MGW_1_2</td><td>0.8750</td><td>0.4405</td></tr><tr><td>MGW_1_3</td><td>0.8435</td><td>0.3410</td></tr><tr><td>MGW_1_4</td><td>0.6547</td><td>0.0510</td></tr><tr><td>N2P_1_1&amp;2</td><td>14.0000</td><td>0.8130</td></tr><tr><td>N2P_1_1&amp;3</td><td>14.4545</td><td>0.8465</td></tr><tr><td>N2P_1_1&amp;4</td><td>12.9091</td><td>0.8395</td></tr><tr><td>N2P_1_2&amp;3</td><td>15.0909</td><td>0.8690</td></tr><tr><td>N2P_1_2&amp;4</td><td>14.0909</td><td>0.8870</td></tr><tr><td>N2P_1_3&amp;4</td><td>12.7273</td><td>0.7815</td></tr><tr><td>H2P_1_1&amp;2</td><td>0.8038</td><td>0.6180</td></tr><tr><td>H2P_1_1&amp;3</td><td>0.8284</td><td>0.7590</td></tr><tr><td>H2P_1_1&amp;4</td><td>0.8024</td><td>0.6750</td></tr><tr><td>H2P_1_2&amp;3</td><td>0.8265</td><td>0.6840</td></tr><tr><td>H2P_1_2&amp;4</td><td>0.8064</td><td>0.5595</td></tr><tr><td>H2P_1_3&amp;4</td><td>0.8246</td><td>0.7040</td></tr><tr><td>V2P_1_1&amp;2</td><td>11.5671</td><td>0.8410</td></tr><tr><td>V2P_1_1&amp;3</td><td>10.4636</td><td>0.7870</td></tr><tr><td>V2P_1_1&amp;4</td><td>10.7841</td><td>0.8240</td></tr><tr><td>V2P_1_2&amp;3</td><td>10.9031</td><td>0.7920</td></tr><tr><td>V2P_1_2&amp;4</td><td>11.6752</td><td>0.8350</td></tr><tr><td>V2P_1_3&amp;4</td><td>8.9590</td><td>0.7120</td></tr><tr><td>FST_1_1&amp;2</td><td>0.0380</td><td>0.6125</td></tr><tr><td>FST_1_1&amp;3</td><td>0.0774</td><td>0.8305</td></tr><tr><td>FST_1_1&amp;4</td><td>0.0603</td><td>0.7725</td></tr><tr><td>FST_1_2&amp;3</td><td>0.0476</td><td>0.8865</td></tr><tr><td>FST_1_2&amp;4</td><td>0.0377</td><td>0.7600</td></tr><tr><td>FST_1_3&amp;4</td><td>0.0428</td><td>0.5855</td></tr><tr><td>LIK_1_1&amp;2</td><td>1.8603</td><td>0.7095</td></tr><tr><td>LIK_1_1&amp;3</td><td>2.2983</td><td>0.8785</td></tr><tr><td>LIK_1_1&amp;4</td><td>2.3294</td><td>0.9120</td></tr><tr><td>LIK_1_2&amp;1</td><td>1.9866</td><td>0.6005</td></tr><tr><td>LIK_1_2&amp;3</td><td>2.1175</td><td>0.7920</td></tr><tr><td>LIK_1_2&amp;4</td><td>2.2919</td><td>0.8830</td></tr><tr><td>LIK_1_3&amp;1</td><td>2.6196</td><td>0.9205</td></tr><tr><td>LIK_1_3&amp;2</td><td>2.3803</td><td>0.9510*</td></tr><tr><td>LIK_1_3&amp;4</td><td>2.3464</td><td>0.8340</td></tr><tr><td>LIK_1_4&amp;1</td><td>2.2651</td><td>0.8720</td></tr><tr><td>LIK_1_4&amp;2</td><td>2.1221</td><td>0.9065</td></tr><tr><td>LIK_1_4&amp;3</td><td>1.9438</td><td>0.5865</td></tr></tbody></table>
<figure><img src="image_40.png" /><figcaption>Figure S1. Identification of the optimal number of genetic clusters (K) inferred by STRUCTURE using the delta(K) and L'(K) methods. (a) and (b) panels correspond to analyses with all 829 individuals, and (c) and (d) panels correspond to STRUCTURE analyses on central Europe dataset.</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>**Figure S2.** Genetic clustering of *C. obsoletus*. (a) Results of the genetic clustering analyses performed with STRUCTURE on complete dataset for K=2 (second most probable K), and (b) Results of the genetic clustering analyses performed with STRUCTURE on dataset of central Europe for K=2 (most probable K). Panels (1) corresponds to an enlargement of the geographical area (a) of the panel circled in red. One specific color has been assigned to each genetic cluster in panels (a) and (b).</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>**Figure S3.** Map of genetic diversity of 2409 cox1 sequence of *C. obsoletus*</figcaption></figure>
<figure><img src="image_43.png" /><figcaption>Figure S4. Neighbor-joining tree of *C. obsoletus* population based on the distance of Cavalli-Sforza & Edwards (1967) computed on microsatellite dataset. The belonging of populations to a genetic cluster defined by the STRUCTURE Analysis is represented by the colors, orange, green and blue for cluster 1, 2 and 3 respectively. Bootstraps values were calculated over 1000 replicates and are represented as percentage and only values superior to 60% are reported.</figcaption></figure>
<figure><img src="image_44.png" /></figure>

**Figure S5.** Heatmap of FST/(1-FST) of all population for microsatellite dataset of *C. obsoletus* in Europe. Levels of genetic distance is represented by color gradient, red mean low level of differentiation and yellow mean high levels of differentiation. The belonging of populations to a genetic cluster defined by the STRUCTURE Analysis is represented by the colors, orange, green and blue for cluster 1, 2 and 3 respectively.
<figure><img src="image_45.png" /><figcaption>Figure S6. PCA summary statistics simulated under scenario 1, the most probable scenarios. Yellow dot represents the real data set. Large dots correspond to simulated posteriors and the small circles correspond to simulated statistics from the prior statistics</figcaption></figure>
References

Abbott RJ, James JK, Milne RI, Gillies AC (2003). Plant introductions, hybridization and gene flow. *Philos Trans R Soc Lond B Biol Sci* **358**(1434): 1123-1132.

Abellan P, Svenning JC (2014). Refugia within refugia – patterns in endemism and genetic divergence are linked to Late Quaternary climate stability in the Iberian Peninsula. *Biol J Linn Soc* **113**: 13–28.

Adamack AT, Gruber B, Dray S (2014). PopGenReport: simplifying basic population genetic analyses in R. *Methods Ecol Evol* **5**(4): 384-387.

Agapow P-M, Burt A (2001). Indices of multilocus linkage disequilibrium. *Mol Ecol Notes* **1**(1-2): 101-102.

Ander M, Troell K, Chirico J (2013). Barcoding of biting midges in the genus *Culicoides*: a tool for species determination. *Med Vet Entomol* **27**(3): 323-331.

Anderson GS, Belton P, Kleider N (1991). *Culicoides obsoletus* (Diptera: Ceratopogonidae) as a causal agent of *Culicoides* hypersensitivity (sweet itch) in British Columbia. *J Med Entomol* **28**(5): 685-693.

Beaumont MA, Rannala B (2004). The Bayesian revolution in genetics. *Nat Rev Genet* **5**(4): 251-261.

Beaumont MA, Zhang W, Balding DJ (2002). Approximate Bayesian computation in population genetics. *Genetics* **162**(4): 2025-2035.

Bertorelle G, Benazzo A, Mona S (2010). ABC as a flexible framework to estimate demography over space and time: some cons, many pros. *Mol Ecol* **19**(13): 2609-2625.

Bhagwat SA, Willis KJ (2008). Species persistence in northerly glacial refugia of Europe: a matter of chance or biogeographical traits? *J Biogeogr* **35**(3): 464-482.

Borkent A (2016). World Species of Biting Midges (Diptera: Ceratopogonidae).

Bourquia M, Garros C, Rakotoarivony I, Gardes L, Huber K, Boukhari I *et al* (2019). Update of the species checklist of *Culicoides* Latreille, 1809 biting midges (Diptera: Ceratopogonidae) of Morocco. *Parasit Vectors* **12**(1): 459.

Bowler DE, Benton TG (2005). Causes and consequences of animal dispersal strategies: relating individual behaviour to spatial dynamics. *Biol Rev Camb Philos Soc* **80**(2): 205-225.

Brito PH, Edwards SV (2009). Multilocus phylogeography and phylogenetics using sequence-based markers. *Genetica* **135**(3): 439-455.

Brookfield JF (1996). A simple new method for estimating null allele frequency from heterozygote deficiency. *Mol Ecol* **5**(3): 453-455.
Calvo JH, Calvete C, Martinez-Royo A, Estrada R, Miranda MA, Borras D *et al* (2009).
Variations in the mitochondrial cytochrome c oxidase subunit I gene indicate northward
expanding populations of *Culicoides imicola* in Spain. *Bull Entomol Res* **99**(6): 583-591.

Carpenter S, Groschup MH, Garros C, Felippe-Bauer ML, Purse BV (2013). *Culicoides* biting midges, arboviruses and public health in Europe. *Antiviral Res* **100**(1): 102-113.

Carpenter S, Wilson A, Mellor PS (2009). *Culicoides* and the emergence of bluetongue virus in northern Europe. *Trends Microbiol* **17**(4): 172-178.

Chakraborty R, Zhong Y, Jin L, Budowle B (1994). Nondetectability of restriction fragments and independence of DNA fragment sizes within and between loci in RFLP typing of DNA. *Am J Hum Genet* **55**(2): 391-401.

Cho HC, Chong CS (1974). Notes On Biting Midges Of The Genus *Culicoides* From South Korea: With Special Reference To Unrecorded Species And Distribution. *Kisaengchunghak Chapchi* **12**(1): 45-75.

Comes HP, Kadereit JW (1998). The effect of Quaternary climatic changes on plant distribution and evolution. *Trends Plant Sci* **3**(11): 432-438.

Cornuet JM, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, Leblois R *et al* (2014). DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. *Bioinformatics* **30**(8): 1187-1189.

Cuellar AC, Kjaer LJ, Kirkeby C, Skovgard H, Nielsen SA, Stockmarr A *et al* (2018). Spatial and temporal variation in the abundance of *Culicoides* biting midges (Diptera: Ceratopogonidae) in nine European countries. *Parasit Vectors* **11**(1): 112.

Delecolle JC (1985). Nouvelle contribution à l'étude systématique et iconographique des espèces du genre *Culicoides*, (Diptéra) : (Cératopogonidae) du Nord-Est de la France. Strasbourg.

Dellicour S, Mardulyn P (2014). spads 1.0: a toolbox to perform spatial analyses on DNA sequence data sets. *Mol Ecol Resour* **14**(3): 647-651.

Despres L (2019). One, two or more species? Mitonuclear discordance and species delimitation. *Mol Ecol* **28**(17): 3845-3847.

Dlugosch KM, Parker IM (2008). Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions. *Mol Ecol* **17**(1): 431-449.

Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR (2014). NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. *Mol Ecol Resour* **14**(1): 209-214.
Earl DA, vonHoldt BM (2011). STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. *Conservation Genetics Resources* **4**(2): 359-361.

Edgar RC (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. *Nucleic Acids Res* **32**(5): 1792-1797.

Ellis JS, Knight ME, Darvill B, Goulson D (2006). Extremely low effective population sizes, genetic structuring and reduced genetic diversity in a threatened bumblebee species, *Bombus sylvarum* (Hymenoptera: Apidae). *Mol Ecol* **15**(14): 4375-4386.

Estoup A, Guillemaud T (2010). Reconstructing routes of invasion using genetic data: why, how and so what? *Mol Ecol* **19**(19): 4113-4130.

Estoup A, Solignac M, Cornuet JM, Goudet J, Scholl A (1996). Genetic differentiation of continental and island populations of *Bombus terrestris* (Hymenoptera: Apidae) in Europe. *Mol Ecol* **5**(1): 19-31.

Evanno G, Regnaut S, Goudet J (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. *Mol Ecol* **14**(8): 2611-2620.

Excoffier L, Smouse PE, Quattro JM (1992). Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. *Genetics* **491**(1985): 479-491.

Frank BM, Piccolo JJ, Baret PV (2011). A review of ecological models for brown trout: towards a new demogenetic model. *Ecology of Freshwater Fish* **20**(2): 167-198.

Fu YX (1997). Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. *Genetics* **147**(2): 915-925.

Funk DJ, Omland KE (2003). Species-Level Paraphyly and Polyphyly: Frequency, Causes, and Consequences, with Insights from Animal Mitochondrial DNA. *Annu Rev Ecol Evol Syst* **34**(1): 397-423.

Gao X, Wang L, Liu T, Xiao J, Wang H (2020). Effect of agro-ecological landscape on the distribution of *Culicoides obsoletus* in northeast China. *Pest Manag Sci*.

Garros C, Balenghien T (2017). Chapitre 14. Les *Culicoides* (Diptera : Ceratopogonidae). In: Éditions I (ed) *Entomologie médicale et vétérinaire*, pp 345-365.

Garza JC, Williamson EG (2001). Detection of reduction in population size using data from microsatellite loci. *Mol Ecol* **10**(2): 305-318.

Goffredo M, Delecolle JC, Semproni G, Lika A (2004). *Culicoides* (Diptera: Ceratopogonidae) in Albania: results of the 2002 entomological survey for bluetongue. *Vet Ital* **40**(3): 209-211.
Gómez A, Lunt DH (2007). *Refugia* within *refugia*: patterns of phylogeographic concordance in the Iberian Peninsula. In: Springer (ed) *Phylogeography of Southern European Refugia*

Goudet J, Raymond M, de Meeus T, Rousset F (1996). Testing differentiation in diploid populations. *Genetics* **144**(4): 1933-1940.

Goulson D, Kaden JC, Lepais O, Lye GC, Darvill B (2011). Population structure, dispersal and colonization history of the garden bumblebee *Bombus hortorum* in the Western Isles of Scotland. *Conserv Genet* **12**(4): 867-879.

Hakima B, Hwang HS, Lee KY (2020). Molecular identification of *Culicoides* (Diptera: Ceratopogonidae) species in Algeria. *Acta Trop* **202**: 105261.

Harrup LE, Laban S, Purse BV, Reddy YK, Reddy YN, Byregowda SM *et al* (2016). DNA barcoding and surveillance sampling strategies for *Culicoides* biting midges (Diptera: Ceratopogonidae) in southern India. *Parasit Vectors* **9**: 461.

Hewitt GM (1999). Post-glacial re-colonization of European biota. *Biol J Linn Soc* **68**(1-2): 87-112.

Hewitt GM (2004). Genetic consequences of climatic oscillations in the Quaternary. *Philos Trans R Soc Lond B Biol Sci* **359**(1442): 183-195; discussion 195.

Hoffmann B, Eschbaumer M, Beer M (2009). Real-time quantitative reverse transcription-PCR assays specifically detecting bluetongue virus serotypes 1, 6, and 8. *J Clin Microbiol* **47**(9): 2992-2994.

Hristescu D, Barbuceanu F, Dascalu L, Nitescu C, Goffredo M, Santilli A *et al* (2020). Species composition and relative abundance of the genus *Culicoides* (Diptera: Ceratopogonidae) in Romania. *Parasit Vectors* **13**(1): 393.

Hubisz MJ, Falush D, Stephens M, Pritchard JK (2009). Inferring weak population structure with the assistance of sample group information. *Mol Ecol Resour* **9**(5): 1322-1332.

Irwin DE, Bensch S, Irwin JH, Price TD (2005). Speciation by distance in a ring species. *Science* **307**(5708): 414-416.

Jacquet S, Garros C, Lombaert E, Walton C, Restrepo J, Allene X *et al* (2015). Colonization of the Mediterranean basin by the vector biting midge species *Culicoides imicola*: an old story. *Mol Ecol* **24**(22): 5707-5725.

Jacquet S, Huber K, Guis H, Setier-Rio ML, Goffredo M, Allene X *et al* (2016a). Spatio-temporal genetic variation of the biting midge vector species *Culicoides imicola* (Ceratopogonidae) Kieffer in France. *Parasit Vectors* **9**: 141.

Jacquet S, Huber K, Pages N, Talavera S, Burgin LE, Carpenter S *et al* (2016b). Range expansion of the Bluetongue vector, *Culicoides imicola*, in continental France likely due to rare wind-transport events. *Sci Rep* **6**: 27247.
James JE, Lanfear R, Eyre-Walker A (2016). Molecular Evolutionary Consequences of Island Colonization. *Genome Biol Evol* **8**(6): 1876-1888.

Kamvar ZN, Tabima JF, Grunwald NJ (2014). Poppr : an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. *PeerJ*.

Kirk H, Dorn S, Mazzi D (2013). Worldwide population genetic structure of the oriental fruit moth (Grapholita molesta), a globally invasive pest. *BMC Ecol* **13**: 12.

Kirkeby C, Bodker R, Stockmarr A, Lind P, Heegaard PM (2013). Quantifying dispersal of european *Culicoides* (Diptera: Ceratopogonidae) vectors between farms using a novel mark-release-recapture technique. *PLoS One* **8**(4): e61269.

Kolleck J, Yang M, Zinner D, Roos C (2013). Genetic diversity in endangered Guizhou snub-nosed monkeys (*Rhinopithecus brelichi*): contrasting results from microsatellite and mitochondrial DNA data. *PLoS One* **8**(8): e73647.

Kosakovsky Pond SL, Frost SD (2005). Not so different after all: a comparison of methods for detecting amino acid sites under selection. *Mol Biol Evol* **22**(5): 1208-1222.

Kozakiewicz CP, Burridge CP, Funk WC, VandeWoude S, Craft ME, Crooks KR *et al* (2018). Pathogens in space: Advancing understanding of pathogen dynamics and disease ecology through landscape genetics. *Evol Appl* **11**(10): 1763-1778.

Lagercrantz U, Ryman N (1990). Genetic Structure of Norway Spruce (Picea Abies): Concordance of Morphological and Allozymic Variation. *Evolution* **44**(1): 38-53.

Lecocq T, Dellicour S, Michez D, Lhomme P, Vanderplanck M, Valterova I *et al* (2013). Scent of a break-up: phylogeography and reproductive trait divergences in the red-tailed bumblebee (*Bombus lapidarius*). *BMC Evol Biol* **13**: 263.

Lecocq T, Rasmont P, Harpke A, Schweiger O (2016). Improving International Trade Regulation by Considering Intraspecific Variation for Invasion Risk Assessment of Commercially Traded Species: The Bombus terrestris Case. *Conservation Letters* **9**(4): 281-289.

Lesieur V, Lombaert E, Guillemaud T, Courtial B, Strong W, Roques A *et al* (2018). The rapid spread of *Leptoglossus occidentalis* in Europe: a bridgehead invasion. *J Pest Sci* **92**(1): 189-200.

Librado P, Rozas J (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. *Bioinformatics* **25**(11): 1451-1452.

Manangwa O, De Meeus T, Grebaut P, Segard A, Byamungu M, Ravel S (2019). Detecting Wahlund effects together with amplification problems: Cryptic species, null alleles and short allele dominance in *Glossina pallidipes* populations from Tanzania. *Mol Ecol Resour* **19**(3): 757-772.

Martinez-de la Puente J, Martinez J, Ferraguti M, Morales-de la Nuez A, Castro N, Figuerola J (2012). Genetic characterization and molecular identification of the bloodmeal sources
of the potential bluetongue vector *Culicoides obsoletus* in the Canary Islands, Spain.
*Parasit Vectors* **5**: 147.

Mathieu B (2011). Les espèces de *Culicoides* du sous-genre Avaritia (Diptera : Ceratopogonidae) dans le monde : révision systématique et taxonomique des espèces d'intérêt dans la transmission d'Orbivirus. 253.

Mayer F, Piel FB, Cassel-Lundhagen A, Kirichenko N, Grumiau L, Okland B *et al* (2015). Comparative multilocus phylogeography of two Palaearctic spruce bark beetles: influence of contrasting ecological strategies on genetic variation. *Mol Ecol* **24**(6): 1292-1310.

Mehlhorn H, Walldorf V, Klimpel S, Jahn B, Jaeger F, Eschweiler J *et al* (2007). First occurrence of *Culicoides obsoletus*-transmitted Bluetongue virus epidemic in Central Europe. *Parasitol Res* **101**(1): 219-228.

Meiswinkel R, van Rijn P, Leijs P, Goffredo M (2007). Potential new *Culicoides* vector of bluetongue virus in northern Europe. *Vet Rec* **161**(16): 564-565.

Mellor PS, Carpenter S, Harrup L, Baylis M, Mertens PP (2008). Bluetongue in Europe and the Mediterranean Basin: history of occurrence prior to 2006. *Prev Vet Med* **87**(1-2): 4-20.

Mignotte A, Garros C, Dellicour S, Gilbert M, Gardès L, Balenghien T *et al* (Submitted in 2020). *Important Dispersal Capacity of One Main Species Vector of Bluetongue and Schmallenberg Viruses Culicoides Obsoletus (Diptera: Ceratopogonidae) Revealed by Landscape Genetic Analyses*.

Mignotte A, Garros C, Gardes L, Balenghien T, Duhayon M, Rakotoarivony I *et al* (2020). The tree that hides the forest: cryptic diversity and phylogenetic relationships in the Palaearctic vector Obsoletus/Scoticus Complex (Diptera: Ceratopogonidae) at the European level. *Parasites & Vectors* **13**(1): 265.

Miranda MA, Borras D, Rincon C, Alemany A (2003). Presence in the Balearic Islands (Spain) of the midges *Culicoides imicola* and *Culicoides obsoletus* group. *Medical and Veterinary Entomology* **17**(1): 52-54.

Murray MD, Kirkland PD (1995). Bluetongue and Douglas virus activity in New South Wales in 1989 : further evidence for long-distance dispersal of the biting midge *Culicoides brevitarsis*. *Austmlian Veterinay Joumal Vol. 72*.

Murrell B, Wertheim JO, Moola S, Weighill T, Scheffler K, Kosakovsky Pond SL (2012). Detecting individual sites subject to episodic diversifying selection. *PLoS Genet* **8**(7): e1002764.

Nei M, Maruyama T, Chakraborty R (1975). The Bottleneck Effect and Genetic Variability in Populations. *Evolution* **29**(1): 1-10.

Nolan DV, Dallas JF, Piertney SB, Mordue Luntz AJ (2008). Incursion and range expansion in the bluetongue vector Culicoides imicola in the Mediterranean basin: a phylogeographic analysis. *Med Vet Entomol* **22**(4): 340-351.
Nychka D, Furrer R, Paige J, Sain S (2017). Fields: Tools for spatial data.

Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin P, O'Hara RB *et al* (2013). *Vegan: Community Ecology Package*. R Package Version. 2.0-10. CRAN.

Onyango MG, Beebe NW, Gopurenko D, Bellis G, Nicholas A, Ogugo M *et al* (2015a). Assessment of population genetic structure in the arbovirus vector midge, *Culicoides brevittarsis* (Diptera: Ceratopogonidae), using multi-locus DNA microsatellites. *Vet Res* **46**: 108.

Pascual M, Chapuis MP, Mestres F, Balanya J, Huey RB, Gilchrist GW *et al* (2007). Introduction history of *Drosophila subobscura* in the New World: a microsatellite-based survey using ABC methods. *Mol Ecol* **16*(15): 3069-3083.

Petit R, Aguinagalde I, de Beaulieu JL, Bittkau C, Brewer S, Cheddadi R *et al* (2003). Glacial refugia: hotspots but not melting pots of genetic diversity. *Science* **300**(5625): 1563-1565.

Phillips BL, Brown GP, Shine R (2010). Life-history evolution in range-shifting populations. *Ecology* **91**(6): 1617-1627.

Pili E, Carcangiu L, Oppo M, Marchi A (2010). Genetic structure and population dynamics of the biting midges *Culicoides obsoletus* and *Culicoides scoticus*: implications for the transmission and maintenance of bluetongue. *Med Vet Entomol* **24**(4): 441-448.

Piry S, Chapuis MP, Gauffre B, Papaix J, Cruaud A, Berthier K (2016). Mapping Averaged Pairwise Information (MAPI): a new exploratory tool to uncover spatial structure. *Methods Ecol Evol* **7**(12): 1463-1475.

Pons O, Petit R (1996). Measuring and testing genetic differentiation with ordered vs. unordered alleles. *Genetics* **144**: 1237-1245.

Pons O, Petit RJ (1995). Estimation, variance and optimal sampling of gene diversity : I. Haploid locus. *Theor Appl Genet* **90**(3-4): 462-470.

Pritchard JK, Stephens M, Donnelly P (2000). Inference of population structure using multilocus genotype data. *Genetics* **155**(2): 945-959.

Provan J, Bennett KD (2008). Phylogeographic insights into cryptic glacial refugia. *Trends Ecol Evol* **23**(10): 564-571.

Prugnolle F, de Meeus T (2002). Inferring sex-biased dispersal from population genetic tools: a review. *Heredity (Edinb)* **88**(3): 161-165.

Prugnolle F, Liu H, de Meeus T, Balloux F (2005). Population genetics of complex life-cycle parasites: an illustration with trematodes. *Int J Parasitol* **35**(3): 255-263.
Pudar D, Petric D, Allene X, Alten B, Ayhan N, Cvetkovikj A *et al* (2018). An update of the *Culicoides* (Diptera: Ceratopogonidae) checklist for the Balkans. *Parasit Vectors* **11**(1): 462.

Purse BV, Carpenter S, Venter GJ, Bellis G, Mullens BA (2015). Bionomics of temperate and tropical *Culicoides* midges: knowledge gaps and consequences for transmission of *Culicoides*-borne viruses. *Annu Rev Entomol* **60**: 373-392.

Purse BV, Mellor PS, Rogers DJ, Samuel AR, Mertens PP, Baylis M (2005). Climate change and the recent emergence of bluetongue in Europe. *Nat Rev Microbiol* **3**(2): 171-181.

Ramos-Onsins SE, Rozas J (2002). Statistical properties of new neutrality tests against population growth. *Mol Biol Evol* **19**(12): 2092-2100.

Rasic G, Schama R, Powell R, Maciel-de Freitas R, Endersby-Harshman NM, Filipovic I *et al* (2015). Contrasting genetic structure between mitochondrial and nuclear markers in the dengue fever mosquito from Rio de Janeiro: implications for vector control. *Evol Appl* **8**(9): 901-915.

Rice WR (1989). Analysing tables of statistic tests. *Evolution* **43**(1): 223-225.

Rodríguez-Sánchez F, Hampe A, Jordano P, Arroyo J (2010). Past tree range dynamics in the Iberian Peninsula inferred through phylogeography and palaeodistribution modelling: A review. *Rev Palaeobot Palynol* **162**(3): 507-521.

Rogers AR, Harpending H (1992). Population growth makes waves in the distribution of pairwise genetic differences. *Mol Biol Evol* **9**(3): 552-569.

Sarvasova A, Goffredo M, Sopoliga I, Savini G, Kocisova A (2014). *Culicoides* midges (diptera: ceratopogonidae) as vectors of orbiviruses in Slovakia. *Vet Ital* **50**(3): 203-212.

Savini G, MacLachlan NJ, Sanchez-Vizcaino JM, Zientara S (2008). Vaccines against bluetongue in Europe. *Comp Immunol Microbiol Infect Dis* **31**(2-3): 101-120.

Schmitt T, Varga Z (2012). Extra-Mediterranean refugia: The rule and not the exception? *Front Zool* **9**(1): 22.

Sghaier S, Hammami S, Goffredo M, Hammami M, Portanti O, Lorusso A *et al* (2017). New species of the genus *Culicoides* (Diptera: Ceratopogonidae) for Tunisia, with detection of Bluetongue viruses in vectors. *Vet Ital* **53**(4): 357-366.

Shao ZY, Mao HX, Fu WJ, Ono M, Wang DS, Bonizzoni M *et al* (2004). Genetic structure of Asian populations of *Bombus ignitus* (Hymenoptera: Apidae). *J Hered* **95**(1): 46-52.

Shults P, Ho A, Martin EM, McGregor BL, Vargo EL (2020). Genetic Diversity of *Culicoides stellifer* (Diptera: Ceratopogonidae) in the Southeastern United States Compared With Sequences From Ontario, Canada. *J Med Entomol* **57**(4): 1324-1327.
Sommer RS, Zachos FE (2009). Fossil evidence and phylogeography of temperate species: 'glacial refugia' and post-glacial recolonization. *J Biogeogr* **36**(11): 2013-2020.

Stewart JR, Lister AM (2001). Cryptic northern refugia and the origins of the modern biota. *Trends Ecol Evol* **16**(11): 608-613.

Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF (1998). Comparative phylogeography and postglacial colonization routes in Europe. *Mol Ecol* **7**(4): 453-464.

Tajima F (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. *Genetics* **123**(3): 585-595.

Takezaki N, Nei M (1996). Genetic distances and reconstruction of phylogenetic trees from microsatellite DNA. *Genetics* **144**(1): 389-399.

Tzedakis PC, Emerson BC, Hewitt GM (2013). Cryptic or mystic? Glacial tree refugia in northern Europe. *Trends Ecol Evol* **28**(12): 696-704.

Wattier R, Engel CR, Saumitou-Laprade P, Valero M (1998). Short allele dominance as a source of heterozygote deficiency at microsatellite loci: experimental evidence at the dinucleotide locus Gv1CT in *Gracilaria gracilis* (Rhodophyta). *Mol Ecol* **7**(11): 1569-1573.

Weir BS, Cockerham CC (1984). Estimating F-Statistics for the Analysis of Population Structure. *Evolution* **38**(6): 1358-1370.

Weir BS, Goudet J (2017). A Unified Characterization of Population Structure and Relatedness. *Genetics* **206**(4): 2085-2103.

White BP, Pilgrim EM, Boykin LM, Stein ED, Mazor RD (2014). Comparison of four species-delimitation methods applied to a DNA barcode data set of insect larvae for use in routine bioassessment. *Freshwater Science* **33**(1): 338-348.

White TA, Perkins SE, Heckel G, Searle JB (2013). Adaptive evolution during an ongoing range expansion: the invasive bank vole (*Myodes glareolus*) in Ireland. *Mol Ecol* **22**(11): 2971-2985.

Whittaker RJ, Field R (2000). Tree species richness modelling: an approach of global applicability? *Oikos* **89**(2): 399-402.

Widmer A, lexer C (2001). Glacial refugia: sanctuaries for allelic richness, but not for gene diversity. *Trends Ecol Evol* **16**(6).

Widmer A, Schmid-Hempel P, Estoup A, Scholl A (1998). Population genetic structure and colonization history of Bombus terrestris s.l. (Hymenoptera : Apidae) from the Canary Islands and Madeira. *Heredity* **81**: 563-572.

Wilson AJ, Mellor PS (2009). Bluetongue in Europe: past, present and future. *Philos Trans R Soc Lond B Biol Sci* **364**(1530): 2669-2681.
Winter DJ (2012). MMOD: an R library for the calculation of population differentiation statistics. *Mol Ecol Resour* **12**(6): 1158-1160.

Yang RC (1998). Estimating Hierarchical F-Statistics. *Evolution* **52**(4): 950-956.
# Chapitre III : Important dispersal capacity of one main species vector of bluetongue and Schmallenberg viruses *Culicoides obsoletus* (Diptera: Ceratopogonidae) revealed by landscape genetic analyses at French scale
A l'échelle continentale, une faible structuration est observée, témoignage d'une histoire évolutive récente combinée à des flux de gènes importants sur l'aire de distribution ouest paléarctique. Mais quels sont les facteurs environnementaux qui conditionnent la connectivité paysagère des *Culicoides* ?

La dispersion des *Culicoides* est connue pour être stratifiée, c'est-à-dire combinant des processus de dispersion actifs, à courte distance, et passifs ou semi-actifs à longue distance, permettant aux individus de faire des bonds de plusieurs centaines de kilomètres. La dispersion active, principalement motivée par la recherche d'hôtes vertébrés pour un repas de sang ou d'un lieu d'oviposition, est conditionnée par les conditions climatiques et météorologiques qui peuvent diminuer l'activité et la survie des *Culicoides*. La dispersion passive par le vent est, elle, influencée par la direction et la vitesse du vent.

Ce troisième chapitre présente pour la première fois une étude de génétique du paysage appliquée aux *Culicoides*. Elle vise à identifier les facteurs environnementaux qui favorisent ou limitent les flux génétiques de *C. obsoletus*. Nos travaux permettent de quantifier l'impact d'un large éventail de facteurs environnementaux et météorologiques classiquement cités comme impactant la dispersion des *Culicoides*. Un plan d'échantillonnage individuel a permis de génotyper 348 individus issus de 46 sites en France à l'aide de 11 marqueurs microsatellites. Une approche de régression multiple sur les matrices de distance nous a permis d'étudier la corrélation entre les distances génétiques et environnementales. L'utilisation d'un raster nul et d'une approche de commonalité nous a permis de nous défaire de la multicolinéarité de ce type d'analyse et d'identifier des contributions uniques et communes des prédicteurs à la variance de la distance environnementale.

Nos travaux permettent de renforcer nos conclusions quant à une dispersion importante des *Culicoides* au-dessus des terres, comme en témoigne l'absence totale d'isolement par la distance à l'échelle française. Une faible différenciation et structure génétique des populations sont observées dans tout le pays. Alors que la contribution unique de tous les facteurs environnementaux et météorologiques testés est faible, celle de la densité du bétail est la plus forte que nous ayons détectée. L'isolement anisotropique significatif sur un axe nord-sud peut être le signe de la dispersion par le vent des *Culicoides*. Une telle capacité de dispersion du *C. obsoletus* sur les territoires terrestres appelle à une ré-évaluation de l'impact des *Culicoides* sur la dispersion du virus. Ce chapitre a été soumis et est sous révision dans le journal *Parasites and Vectors*. https://dx.doi.org/10.21203/2Frs.3.rs-75603%2Fv1
# Important dispersal capacity of one main species vector of bluetongue and Schmallenberg viruses *Culicoides obsoletus* (Diptera: Ceratopogonidae) revealed by landscape genetic analyses at French scale

Antoine Mignotte<fnref n="1" />,<fnref n="2" />*, Claire Garros<fnref n="1" />,<fnref n="2" />, Simon Dellicour<fnref n="3" />,<fnref n="4" />, Maude Jacquot<fnref n="3" />,<fnref n="5" />, Marius Gilbert<fnref n="3" />, Laetitia Gardès<fnref n="1" />,<fnref n="6" />, Thomas Balenghien<fnref n="1" />,<fnref n="7" />,<fnref n="8" />, Maxime Duhayon<fnref n="1" />,<fnref n="2" />, Ignace Rakotoarivony<fnref n="1" />,<fnref n="2" />, Maïa de Wavrechin<fnref n="1" />,<fnref n="2" />, Karine Huber<fnref n="1" />.

### Emails:

AM: antoine.mignotte@cirad.fr

CG: claire.garros@cirad.fr

SD: sdellico@ulb.ac.be

MJ: maude.jacquot@inrae.fr

MG: marius.gilbert@ulb.ac.be

LG: laetitia.gardes@cirad.fr

TB: thomas.balenghien@cirad.fr

MD: maxime.duhayon@cirad.fr

IR: ignace.rakotoarivony@cirad.fr

MW: maiadewav@gmail.com

KH: karine.huber@cirad.fr

<footnote n="1">¹ ASTRE, Univ Montpellier, Cirad, INRAE, Montpellier, France.</footnote>

<footnote n="2">² Cirad, UMR ASTRE, F-34398 Montpellier, France.</footnote>

<footnote n="3">³ Spatial Epidemiology Lab (SpELL), Université Libre de Bruxelles, CP160/12, 50, av. FD Roosevelt, 1050 Bruxelles, Belgium</footnote>

<footnote n="4">⁴ Department of Microbiology, Immunology and Transplantation, Rega Institute, KU Leuven, Herestraat 49, 3000 Leuven, Belgium</footnote>

<footnote n="5">⁵ Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA, F-63122 Saint-Genès-Champanelle, France</footnote>

<footnote n="6">⁶ Cirad, UMR ASTRE, F-97170 Petit-Bourg, Guadeloupe, France.</footnote>

<footnote n="7">⁷ Cirad, UMR ASTRE, M-10100 Rabat, Morocco.</footnote>

<footnote n="8">⁸ Institut Agronomique et Vétérinaire Hassan II, Unité Microbiologie, immunologie et maladies contagieuses, 10100 Rabat-Instituts, Morocco.</footnote>

<footnote>* Correspondence: antoine.mignotte@cirad.fr, claire.garros@cirad.fr and karine.huber@cirad.fr</footnote>
**Background:** In the last two decades, recurrent epizootics of bluetongue virus and Schmallenberg virus have been reported in the western Palearctic region. These viruses affect domestic cattle, sheep, goats and wild ruminants and are transmitted by native hematophagous midges of the genus *Culicoides* (Diptera: Ceratopogonidae). *Culicoides* dispersal is known to be stratified, *i.e.* due to a combination of dispersal processes occurring actively at short distances and passively or semi-actively at long distances, allowing individuals to jump hundreds of kilometers.

**Methods:** Here, we aim to identify the environmental factors that promote or limit gene flows of *Culicoides obsoletus*, an abundant and widespread vector species in Europe using an innovative framework integrating spatial, population genetics and statistical approaches. A total of 348 individuals were sampled in 46 sites in France and have been genotyped using 11 newly-designed microsatellite markers.

**Results:** We found a low genetic differentiation and population structure for *C. obsoletus* across the country. Using three complementary inter-individual genetic distances, we did not detect any significant isolation-by-distance, but a significant anisotropic isolation-by-distance on a North-South axis. We have employed a multiple regression approach on distance matrices to investigate the correlation between genetic and environmental distances. Among all the environmental factors that have been tested, only cattle density seems to have an impact on *C. obsoletus* gene flow.

**Conclusions:** Such a high dispersal capacity of *C. obsoletus* over land territories calls for a re-evaluation of the impact of *Culicoides* on virus dispersal and highlights the urgent need to better integrate molecular, spatial and statistical information to guide vector-borne diseases control.

**Keywords:** *Culicoides obsoletus*, landscape genetics, microsatellites, dispersal, Palearctic region
## Background

During a vector-borne disease outbreak, a precise understanding of the dispersal capacity of the vector species is a key aspect to implement appropriate control strategies in order to limit the spread of the disease. Dipteran dispersal is known to be impacted by the composition and configuration of the landscape. When the dispersal of an organism is favored by a landscape characteristic (“environmental factor”) it can be categorized as a conductance factor (Manel and Holderegger, 2013). On the contrary, when the dispersal is limited by a landscape characteristic, it can be categorized as a resistance factor (Manel et al, 2007; Manel and Holderegger, 2013). However, some types of dispersal allow overcoming resistance factor, as it is the case for “stratified dispersion”. *Culicoides* dispersal is described as “stratified”, due to a combination of dispersal processes occurring actively at short distances, and passively or semi actively at long distances (Murray and Kirkland, 1995).

Northern Europe experienced very sudden and rapid outbreaks of bluetongue virus (BTV) in 2006-2008 and of Schmallenberg virus (SBV) in 2011-2012 (Carpenter et al, 2013). Both viruses have spread very quickly and massively across the whole western Palearctic region transmitted by native *Culicoides* species. *Culicoides obsoletus* has been identified as the main vector species responsible for the transmission of BTV and SBV to wild and domestic ruminants in western Palearctic region (Carpenter et al, 2009; Hoffmann et al, 2009; Savini et al, 2008). *Culicoides obsoletus* is a very widespread species able to breed in a wide range of habitats (Mehlhorn et al, 2007; Purse et al, 2015). BTV spread is facilitated by favorable conditions for midge activity and viral replication and the environmental and climatic drivers of BTV transmission in Europe have been suggested (Cuellar et al, 2020; Guis et al, 2012; Jacquot et al, 2017; Maclachlan and Guthrie, 2010; Purse et al, 2005). For example, high ambient temperature reduce incubation period of the virus (Wittmann et al, 2002) whereas high precipitation and wind speeds can reduced the *Culicoides* flight activity (Baylis et al, 2010; Carpenter et al, 2008c; Gloster et al, 2008).

Previous studies of mark-release-recapture on *Culicoides* species showed that the post-release dispersal distance travelled, during two nights, ranges from 1 to 2.5 km, and is linked to the gradual search for hosts or oviposition sites (Kluiters et al, 2015b). The maximum recapture travelled distance recorded was 6 km for *Culicoides mohave* in a particular desert landscape (Brenner et al, 1984). Yet, the potentially high dispersal capacity of *Culicoides* and high mortality of adults when manipulating represents a practical limitation in the context of a mark-
release-recapture procedure. Dispersal phenomena driven by wind currents on a much larger
scale have also been detected during outbreaks related to *Culicoides-borne diseases* (Alba *et al*, 2004; Ducheyne *et al*, 2007; Eagles *et al*, 2012; Gloster *et al*, 2008; Hendrickx *et al*, 2008; Jacquet *et al*, 2016b; Murray, 1987; Murray and Kirkland, 1995). Introduction of BTV
serotypes by wind-borne infected *Culicoides* has been demonstrated from Northern Africa to
Southern Europe (Gubbins *et al*, 2008), from Kenya to South-Western Indian Ocean Islands
(Onyango *et al*, 2015a), from Sardinia to the Balearic Islands (Alba *et al*, 2004), from France
to Corsica (Jacquet *et al*, 2016a), in Ireland (McGrath *et al*, 2018), or in UK (Burgin *et al*,
2017). Most of these studies are supported by modelling analyses in which dispersal trajectories
have been evaluated by atmospheric dispersion models over water bodies (Agren *et al*, 2010;
Aguilar-Vega *et al*, 2019; Andrew Jones, 2007; Durr *et al*, 2017; Eagles *et al*, 2014; McGrath
*et al*, 2018). *Culicoides* dispersal over land territories has been less investigated (Burgin *et al*,
2013) and little attempt has been made to link it to environmental factors. Given this high
dispersal capacity and the stratified dispersal pattern of *Culicoides*, it is crucial to determine
inland connectivity among populations and to identify the potential environmental factors that
promote or limit gene flows between them.

Classic approaches to study the impact of landscape on gene flow generally use
population-based sampling and are often based on a limited number of sampled populations
(Jacquet et al, 2015; Onyango et al, 2015b). However, *Culicoides* such as other dipterans
species, are not spatially structured into separate populations and must be considered as a
continuum of individuals heterogeneously distributed across landscape. In order to identify the
environmental factors having an impact on the dispersal of *Culicoides*, an individual approach
is more relevant as it avoids inter-population genetic distance misinterpretations (Shirk et al,
2017). Indeed, the individual approach in landscape genetics aims at maximizing the number
of sampling sites, and thus brings a much greater statistical power to detect spatial patterns of
genetic differentiation and the environmental factors that cause them (Landguth et al, 2010). In
addition to the analysis of inter-individual genetic distances, population genetic structure can
also be investigated, e.g. using gold-standard Bayesian clustering methods. However, if
isolation-by-distance occurs, these clustering analyses could overestimate the actual number of
genetic clusters (Frantz et al, 2009). It is therefore important to also incorporate a visualization
approach into the analytical workflow (Piry et al, 2016). MAPI allow to visually compare
genetic dissimilarity with some environmental factors but also to develop working hypotheses.
On the other side, it is also necessary to use as a complement statistical analyses of genetic and
environmental distances.

The present study here aims at determining population inland connectivity at large
geographical scales of a main vector species and to identify the environmental factors that
promote or limit gene flows between *Culicoides* populations in France. For this purpose, we
characterized 13 microsatellite markers dedicated to *C. obsoletus*. We are proposing a
complementary framework integrating multiple approaches which can be applied more
generally to the study of gene flow and links with environmental factors.
# Methods

## *Culicoides obsoletus* sampling and morphological identification

The data set analysed in the present study is composed of 368 individual females biting midges collected in 46 sites in France in April 2011, using the national surveillance network for *Culicoides* populations or local complementary collections (Table S1). Collections were made overnight using Onderstepoort Veterinary Institute (OVI) light traps set up at farms near cattle or sheep. All insects were stored in 70% ethanol. Morphological identification to species level has been carried out in the same way as described in a previous study (Mignotte *et al*, 2020) under a binocular microscope using the available identification keys (Delecolle, 1985).

## DNA extraction, amplification, sequencing and sequence analyses

A total of 368 individuals belonging to the *C. obsoletus* species complex were DNA-extracted, cox1 amplified and sequenced following the procedure detailed in Mignotte *et al*. (2020). To ensure that only *C. obsoletus* individuals were included in the final data set (Garros *et al*, 2014), sequences assignation was run according to Mignotte *et al*. (2020). The reference cox1 sequences used to make a specific assignment of all individuals are available in Table S2. All cox1 sequences were then aligned using the MUSCLE algorithm (Edgar, 2004) implemented in the software package GENEIOUS v.6.0.5 (Biomatters, http://www.geneious.com). In order to estimate unbiased gene flows of *C. obsoletus*, only individuals identified as such were used for further analyses.

## Development of microsatellite markers

Microsatellite markers were defined by next-generation sequencing following a similar protocol to the one described by Morillo *et al*. (2016). *Culicoides obsoletus* samples from France and United Kingdom (Table S3) were used as biological material for the DNA library. Among the 82 276 simple sequence repeats (SSRs) identified, primers were designed for 43 676 SSRs, of which 15 888 were dinucleotides, 23 656 were trinucleotides, and 4,132 contained complex SSR motifs. For the SSR screening, we selected 3 135 primer pairs flanking dinucleotide SSR motifs with a minimum of 12 repetitions and amplifying fragments between 150 and 300 base pairs (bp) in length and a half denaturation temperature (Tm) close to 55 °C. We screened a set of 30 primer pairs and optimized both polymerase chain reaction (PCR) and multiplex PCR conditions.
Microsatellite genotyping and loci filtering

Microsatellite makers were amplified by multiplex PCR with the Type-it-Microsatellites kit (Qiagen, Valencia, CA, USA) according to the protocol described in the manufacturer's manual and the annealing temperature given in Table S4. Standard conditions for PCR amplification included an initial denaturation step at 95 °C for 5 min, 35 cycles of denaturation for 30s at 95 °C, annealing for 1 min at variable temperature (Table S4) and elongation for 1 min at 72 °C, followed by a final elongation of 5 min at 72°C. Fragments were separated on an Applied Biosystems 3500xL Genetic Analyzer. Allelic size allocations for all individuals and microsatellite markers were performed using the program GeneMapper® version 5 (Life Technologies, Applied Biosystem) with double blind reading to limit the potential interpretation bias from the reader.

A first validation of the polymorphism of these 30 markers was carried out by amplifying and genotyping 48 individuals from various locations within the species geographic range. Then 13 loci showing good genetic profiles and clear allelic size variability were characterized as polymorphic markers and were selected for further analyses (Table S4).

The 368 *C. obsoletus* previously extracted were then amplified and genotyped at these 13 microsatellite loci. Genetic variability parameters, such as the number of alleles per locus, the allele size range, the observed and expected heterozygosity were estimated per locus, and per population over the entire dataset using the R packages *ggplot* (Ginestet, 2011), *poppr* (Kamvar *et al*, 2014) and *mmod* (Winter, 2012). Short allele dominance could be a source of heterozygote deficiency in a microsatellites data set (Manangwa *et al*, 2019; Wattier *et al*, 1998). In order to avoid such bias, Spearman's correlations rank between the size of the alleles and their frequencies were calculated in R (Wattier *et al*, 1998). Null allele frequencies were estimated using the R package *PopGenReport* (Adamack *et al*, 2014), and a Bonferroni's correction was applied to all matched tests to take into account potential biases associated with multiple comparisons (Rice, 1989). Correlation between *F*<sub>IS</sub> and null allele was calculated in R, and locus with significative correlation were removed. The presence of an imbalance in genetic binding between each pair of loci has been tested by accurate Fisher tests performed with the R package *poppr* (Kamvar *et al*, 2014). The same R package was used along with the *pegas* R package to identify a potential gap in the panmictic matching regime. For this purpose, we performed a chi-square tests to compare the observed heterozygosity deficits (Ho) to the expected heterozygosity (He) under the Hardy-Weinberg (HW) equilibrium with a significance threshold of 5%. We also calculated the fixation index (*F*<sub>IS</sub>) (Weir and Cockerham, 1984).
## Genetic clustering analysis

In order to determine the most likely number of genetic clusters (K) and the probability of assignment of each individual to these clusters, we used the Bayesian clustering methods implemented in the program STRUCTURE (version 2.3.4) (Pritchard *et al*, 2000) and GENELAND (version 4.5.0) (Guillot *et al*, 2005). The algorithm implemented in STRUCTURE infers genetic clusters that minimize deviations from Hardy-Weinberg equilibrium. In STRUCTURE, we specified the following settings: correlated allele frequency, an admixture model, and a “locprior” model. These options respectively allow to: assume allele frequencies are similar among populations, estimate for each individual the belonging proportion to each K clusters, and use sampling locations as prior information to assist the clustering, which is particularly adapted when there is a weak genetic structuring (Hubisz *et al*, 2009). We performed 10 independent runs for each value of K varying from 1 to 10. Each run consisted in a burn-in of 100,000 iterations in the Markov Monte Carlo chain (MCMC) and a stationary phase of 1,000,000 iterations in the MCMC. The most probable number of clusters was inferred with the ΔK method (Evanno et al. 2005), via the STRUCTURE Harvester web platform (Earl and vonHoldt, 2011). This ΔK value is a way to determine the inflection in the [Ln P(D)] curve. The clustering resulting from the Bayesian inference was transposed into percentages of assignment of each individual to the K inferred clusters and plotted on a map.

GENELAND also uses a Bayesian algorithm to infer population genetic clusters while taking into account the spatial position of individuals, making it a spatially-explicit clustering method. The most probable number of clusters was also determined by running the algorithm with K ranging from 1 to 10. The analysis was based on 1,000,000 MCMC iterations with a thinning of 1000, maximum rate of the Poisson process fixed to 100, maximum number of nuclei in the Poisson-Voronoi tessellation fixed to 300 and a burn-in of 100. We used the R package *graphics* to produce a distribution map of genetic structure resulting of STRUCTURE and GENELAND analyses.

## Computation of genetic distances and mapping genetic dissimilarity

We computed three complementary inter-individual genetic distances. The first inter-individual genetic distance is based on a Factorial Correspondence Analysis (FCA) performed with the program GENETIX v. 4.05.2 (Belkhir, 2004). The first ten FCA axes were used to calculate an Euclidean distance between all individuals using the R package *Ecodist* (Goslee and Urban,
2007), hereafter called the “FCA” distance (Kimmig *et al*, 2020). The second inter-individual
genetic distance is the Rousset’s distance (a<sub>R</sub>), (Rousset, 2000), corresponding to the F<sub>ST</sub>/(1
–F<sub>ST</sub>) but for pairs of individuals. We used the SPAGeDi 1.5 program (Hardy and Vekemans,
2002) to compute a<sub>R</sub> genetic distance. The third inter-individual genetic metric is the Loiselle’s
kinship coefficient (“LKC”, (Loiselle *et al*, 1995)), which, unlike FCA and a<sub>R</sub> that are measures
of genetic dissimilarity, is an index of similarity (Hardy, 2003). We used the R package *fields*
(Nychka *et al*, 2017) to investigate isolation-by-distance by performing linear regressions and
Mantel tests between inter-individual genetic distances and log-transformed geographic
distances with R package *fields* (Nychka *et al*, 2017). Mantel tests (Mantel, 1967) were
performed with the R package *vegan* (Oksanen *et al*, 2013). Inter-individual genetic distances
were mapped using the program MAPI (Piry *et al*, 2016). The method implemented in MAPI
allows to visualize spatial variation in genetic dissimilarity. It can also define areas where
genetic dissimilarity is significantly lower or higher than expected by chance (Piry *et al*, 2016).
MAPI analyses were based on 1,000 permutations and an α value of 0.05.

Impact of environmental factors on genetic differentiation

We investigated the potential impact of several environmental factors on the genetic differentiation of *C. obsoletus*. In practice, we tested the association between matrices of pairwise genetic distances (see above) and matrices of environmental distances. These environmental distances were computed using the program CIRCUITSCAPE 4.0.5 (McRae and Shah, 2008; McRae, 2006) and were based on distinct environmental rasters (i.e. geo-referenced grids of environmental values) either treated as potential resistance (R) or conductance (C) factors:

i) host densities: European index of distribution of roe deer (C) (Neil *et al*, 2014), European index of distribution of red deer (C) (Neil *et al*, 2014), global cattle, sheep and goat distributions in 2010 (C) (Gilbert *et al*, 2018).

ii) terrestrial habitat: mean global elevation in 2010 (R) (Carabajal *et al*, 2011) and type of landscape cover (urban area, grassland, forest areas and croplands) (R). The four distinct land cover rasters (resolution: ~1000 m) were generated from the Corine Land Cover 2012 raster (CLC; www.eea.europa.eu; resolution: ~100 m) (Fig. S1).

iii) meteorological and climate data: mean surface temperature (R), mean precipitation (R), and mean wind speed and direction (C). We used two
alternative rasters for each variable: April 2011 mean and 2000-2010 mean. Surface temperatures were obtained from monthly day and night land surface temperature mean (Wan, 2014). Precipitation and wind rasters were obtained from monthly means of precipitation, wind speed and the U and V components of wind (components of the horizontal wind towards the east and North) from the European Centre for Medium-range Weather Forecasts (Era5 Interim, http://www.ecmwf.int) over 10 years (2001-2010). Wind U and V components were used to compute prevailing meteorological wind direction (in degrees). Map and arrowheads indicating the mean wind direction over the 10 years at each raster pixel were generated with the R package *Rwind* (Fig. S2) (Fernández-López and Schliep, 2018).

A so-called “null raster” was also created and used as a negative control, i.e. a raster on which environmental distances computations correspond to a proxy of geographical distances. The null raster is a uniform raster with all cell values equal to 1. As in Dellicour et al. (2019), three values of *k* (10, 100 and 1000) were used to modify the potential impact of the environment on the resistance/conductance value. We used multiple regressions on distance matrices (MRDM) coupled with commonality analyses (CA; (Newton and Spurrell, 1967)) to identify unique and common contributions of predictors to the variance in the environmental distance (response variable). All the analyses were carried out with the R packages *ecodist* (Goslee and Urban, 2007) and *yhat* (Nimon et al., 2008). In addition, we also performed univariate analyses by comparing (i) the determination coefficient $R^2$ obtained from the linear regression of genetic distances on distances computed on the environmental raster and (ii) the determination coefficient $R^2$ obtained from the linear regression of genetic distances on distances computed on the null raster. Only the environmental factors associated with a $R^2$ higher than that the one obtained from the linear regression based on environmental distances computed on the null raster were selected for the multivariate analyses (MRDM-CA) (Dellicour et al., 2019). The same criterions to identify a suppressors were used for the multivariate analysis as in the commonality analyses descriptions [79, 80]. These suppressors allow to select environmental variables that are more explanatory than the geographical distance alone. An environmental variable is considered to be a suppressor if regression coefficient and correlation, or when the unique contribution and the common contribution are of opposite sign (Prunier et al., 2017).
## Anisotropic isolation-by-distance

In order to identify directional gene flows, we must calculate the projected distance between two sampling sites. The angle between sites was calculated using the R package geosphere (Hijmans, 2019) and transformed from degrees to radians. The set of angles between 0 and 360 were then tested as angles that maximize isolation-by-distance. The great-circle geographic distance between sampling sites was computed with the R package fields and then log-transformed (Nychka et al, 2017). The projected distance matrix was calculated for each angle between all sampling sites, using the formula below. In calculation of the projected geographical distance as a function of the angle between populations, dAB is the geographical distance between population A and B and aAB is the angle between populations A and B.

$$log ( | d_{AB} * cos(a_{AB}) | )$$

A linear regression of the genetic distances on the projected geographical distances obtained for each angle was then performed. The angle that will maximize the R² of this regression with a positive regression coefficient was considered as the angle maximizing the isolation-by-distance signal (Jay et al, 2013). In order to confirm this analysis with a similar but complementary approach, we evaluated the direction that maximizes the correlation between geographic and genetic distance using the PASSAGE program and a "Bearing" analysis. The workflow summarizing the genetic analyses carried-out is described in Figure S3.
# Results

## Characterization of microsatellite markers and genetic diversity

A total of 368 cox1 sequences of *C. obsoletus* were obtained after the DNA barcoding step (Table S3) and 368 *C. obsoletus* from 46 populations were genotyped with 13 microsatellite loci (Table S4). A filter has also been applied to exclude missing data. Individuals with more than 10% missing data were excluded from the analysis, i.e. 20 individuals. Loci OBSms25 and OBSms26 had the highest proportion of null alleles with 56% and 22% respectively, but also the highest FIS values of 0.691 and 0.357 respectively (Table S5). The expected heterozygosity, which varies from 0.918 (OBSms5) to 0.404 (OBSms29) is always higher than that observed, varying from 0.872 (OBSms13) to 0.197 (OBSms25), reflecting a heterozygous deficit (Table S5). To ensure that all markers provide independent information, linkage disequilibrium between each pair of loci was tested. After Bonferroni correction, no significant linkage disequilibrium was observed between all loci. The OBSms25 locus was also the only locus with more than 5% missing data. Thus, only markers with a null percentage inferior to 10% were kept for the final analysis. In view of the high proportion of null alleles and the strong correlation between FIS and the rate of null alleles (r=0.947, p-value>0.001), the OBSms25 and OBSms26 loci were excluded from the analysis. No significant correlation between allele size and allele frequency was found for any of the 11 remaining loci (Table S5). The genetic resolution of microsatellite markers depends on the number of markers and their polymorphism. All selected 11 loci were therefore polymorphic (Table S5) and the number of alleles present at a given locus varied from 40 for OBSms5 to 9 for OBSms29. No significant isolation-by-distance was detected with the different inter-individual genetic distances, aR (R²<0.001, p=0.439, Fig. S4), LKC (R²<0.001, p=0.992, Fig. S4), and FCA (R²<0.001, p=0.361, Fig. S4).

## Clustering analyses and mapping of genetic distances

The Bayesian analysis performed with STRUCTURE revealed a low genetic structuring of *C. obsoletus*. Meanwhile, the exact number of genetic clusters was not clearly defined, the optimal value of K was 2 (Fig. S5); i.e. the co-existence of two genetic clusters was the most likely statistically (Fig. 1). The assignment of individuals to the different clusters revealed a weak genetic structuring of the populations. Although the assignment of a total of 284 out of the 348 *Culicoides* individuals to a genetic cluster is well supported (Q>0.90), no clear association between spatial clustering and distribution can be observed, i.e. there is no association between
*Culicoides* sampling location and the genetic clusters distributions. However, a West-East transect seems to emerge as samples collected on this transect were more frequently assigned to cluster 1, contrary to the North-South transect more represented by cluster 2. Genetic clusters inferred by GENELAND and STRUCTURE (K=2) are consistent with these results, both clusters being found at all the sampling sites (Fig. 1 ab). Results for the second optimal K value (K=4) of the STRUCTURE program gave a similar view, reinforcing the conclusion that populations are genetically weakly structured (Fig. S6). Additional STRUCTURE analysis has been carried out with only loci OBSms 4, 11 and 28 which are at Hardy Weinberg's equilibrium and no additional structure is detected using only this loci (data not shown). Genetic dissimilarity result from MAPI analyses with the $a_R$ genetic distance showed a non-homogeneous pattern of dissimilarity. The maximum inter-individual genetic dissimilarity was of 0.182 and the lowest was of 0.002. The South of France was the most genetically dissimilar region. The center of France was also a zone of strong inter-individual dissimilarity, while eastern France showed low genetic dissimilarity. Another area spanning from the Paris region to Alsace (North east of France) showed low genetic dissimilarity (Fig. 1 cd). MAPI analyses with FCA inter-individual genetic distance, showed a homogeneous genetic dissimilarity. However, almost none of the areas was significantly genetically dissimilar in MAPI analysis with the three genetic distances used.

Investigating the impact of environmental factors

Univariate regressions revealed a significant positive association of FCA genetic distances with cattle distributions for *k*=100 (Q=0.009, p-value=0.001) and *k*=1000 (Q=0.029, p-value=0.001), with sheep distributions for *k*=100 (Q=0.002, p-value=0.002) and for *k*=1000 (Q=0.006, p-value=0.001), with elevation for *k*=100 (Q=0.001, p-value=0.005) and for *k*=1000 (Q=0.001, p-value=0.006), and with grassland for *k*=100 (Q=0.009, p-value=0.001) and for *k*=1000 (Q=0.011, p-value=0.001) distances (Table S6). Similarly, $a_R$ inter-individual genetic distances were significantly associated with distances computed from elevation for *k*=100 (Q=0.00018, p-value=0.006) and for *k*=1000 (Q=0.00038, p-value=0.007), with grassland for *k*=100 (Q=0.001, p-value=0.001) and for *k*=1000 (Q=0.002, p-value=0.001) and with urban areas for *k*=100 (Q=0.004, p-value=0.001) and for *k*=1000 (Q=0.007, p-value=0.001) rasters (Table S6). However, no significant association was found between “LKC” genetic distance and the environmental distances (Q values always <0).
Multivariate regressions revealed a significant positive association of FCA inter-individual genetic distance with cattle distributions for *k*=1000 (*r*=0.1733, p-value <0.001) and with grassland for *k*=1000 (*r*=0.1080, p-value <0.001) (Table S6). Also, multivariate regressions did not reveal any significant positive association of *a*<sub>R</sub> inter-individual genetic distances with environmental distances. All environmental factors were considered as suppressors when using *a*<sub>R</sub> to measure inter-individual genetic distances. Cattle distribution (*k*=1000) with FCA inter-individual genetic distance was the only environmental factor that showed a unique contribution to the variance in the dependent variable higher than 1% (U=0.02) (Table S6).

## Anisotropic isolation-by-distance

Globally, correlation between the *genetic distance* and the distance along the bearing θ, *d*<sub>θ</sub>, changed as a function of bearing θ. Anisotropic isolation-by-distance analyses revealed a significant positive correlation between *a*<sub>R</sub> and FCA inter-individual genetic distances and orientational distances computed along North–South orientation and a significant negative correlation between genetic distances and orientational distances computed along a West-East orientation (Fig. 2 a,b). For the analysis with LKC, which is a variable of genetic similarity, the North–South was the orientation that maximized the isolation-by-distance (Fig. 2 c). Bearing analyses with PASSAGE software revealed a significant distance isolation on the North/South axis with the genetic distances *a*<sub>R</sub> and FCA (Fig. 3 a,b). A negative correlation between LKC genetic similarity and geographical distance was observed on the North/South axis (Fig. 3 c).
Discussion

This study presented for the first time an extensive landscape genetic study carried out on a
main vector species of interest for livestock diseases, *C. obsoletus*. Bayesian clustering analysis
revealed a very weak genetic structure of *C. obsoletus* in France and a low level of inter-
individual genetic differentiation was observed. We assessed the impact of a wide range of
environmental factors on this pattern. Univariate and multivariate analyses highlighted an
absence, or a weak impact of most of the tested environmental factors on inter-individual
measures of pairwise genetic differentiation. No isolation-by-distance pattern was detected at
the French scale with the set of inter-individual genetic distances used. However, an anisotropic
isolation-by-distance analysis revealed a significant distance isolation on a North-South axis.

Important gene flows between *C. obsoletus* populations

The low overall inter-individual genetic dissimilarity highlighted by our study reflected the
important level of gene flow for *C. obsoletus* and comes to reinforce what has been described
previously for other *Culicoides* species. In particular, important level of gene flows have
already been observed in France for *C. imicola* (Jacquet *et al*, 2016a) and in Australia for *C.
brevitarsis* (Onyango *et al*, 2015a). Genetic studies on *C. imicola* in Europe revealed an
important gene flow as reflected by the inference of two large genetic clusters: a “Central
Mediterranean cluster” including Algeria, Sardinia, Corsica, Pyrénées-Orientales and Var
French departments of France, and a “Western Mediterranean cluster” including Morocco,
Spain, Portugal and Majorca (Jacquet *et al*, 2016b). In North America, similar results were also
found in *C. stellifer*. No barriers to gene flow could be identified in this species in the southeast
United States (Shults *et al*, 2020). In addition, no isolation-by-distance has been observed at the
French scale. Identifying a genetic structure pattern is a real challenge when inter-individual
genetic dissimilarity is low. Although STRUCTURE performs well at low levels of population
differentiation (0.02 < FST < 0.10) by using prior population and correlated allele frequency
models (Latch *et al*, 2006), when the differentiation is weaker, such as in the case of highly
dispersed organisms, the latter may encounter difficulties (Duchesne and Turgeon, 2012). This
underlines the importance of using complementary tools of visualization approaches like
MAPI. In addition, the use of an individual approach avoids to artificially consider individuals
grouped into populations solely due to the effect of the trapping sites (Manel *et al*, 2007; Prunier
*et al*, 2013a). Although our method of analysis seems complementary and coherent to detect a
genetic structure, *C. obsoletus* is not genetically or geographically structured at the French
scale.
Livestock densities as potential drivers of *C. obsoletus* gene flow

It is undeniable that large gene flows homogenize the genetic diversity of *Culicoides*. This can be explained by active dispersal and host-seeking movements, as suggested by our results. While the unique contribution of all environmental factors tested was very small, the one of cattle density is the strongest we have detected. Host density as a conductance factor for *Culicoides* is obviously a factor consistent with the biology of the species. These results are in line with those obtained in landscape genetics studies of the BTV, which had identified distributions of cattle and sheep as a key factors in BTV dispersal (Jacquot *et al*, 2017). In addition, previous studies showed that dairy cattle density is negatively correlated with BTV spread. Although paradoxical with the previous conclusions, this could be explained by the fact that dairy cattle are clustered around the milking parlour and move little, forming a fixed feed source, limiting the spread of *Culicoides*. On the contrary, beef cattle disperse much more in the pastures and by this phenomenon could encourage the diffusion of *Culicoides* in search of a source of blood (Pioz *et al*, 2012). In view of these converging results, of the marked preference of certain species of *Culicoides* vector for cattle (Viennet *et al*, 2013), and of BTV emergence or reemergence events in cattle in the Netherlands in 2006 and in France in 2015 (Elbers *et al*, 2008; Sailleau *et al*, 2017), it seems crucial to closely monitor the surveillance of *Culicoides* in the vicinity of beef cattle farms (Viennet *et al*, 2013).

In addition, it has been shown that BTV spread is facilitated at low elevation, up to 300 m (Jacquet *et al*, 2015; Pioz *et al*, 2012). Altitude, an environmental factor tested in this study, failed to explain the inter-individual genetic differentiation of *C. obsoletus*. The altitude therefore does not act as a barrier between sampling sites located at relatively low altitude. It would thus be interesting to integrate in future studies sampling sites with a high altitude above 1,000 m, which is possible for *C. obsoletus* highlighting a large altitudinal range (Diaz-Sanchez *et al*, 2018). However, it is possible that the altitude does not directly impact the dispersal activity of *Culicoides*, but only the replication or viral infection of BTV due to low temperature. *Culicoides obsoletus* is extremely generalist and can take a blood meal on a wide range of hosts (Elbers and Meiswinkel, 2014). The study of phylogenetically close but ecologically very different species of *Culicoides*, such as *C. chiopterus* (Ayllon *et al*, 2014), might then allow the identification of very different dispersal patterns and an importance of bovine density, the latter being exclusively dependent on it for egg-laying (Dijkstra *et al*, 2008).
## Long-dispersal of *C. obsoletus*

The most genetically dissimilar individuals are mainly located in the southernmost populations of the sampling area. Multiple non-exclusive lines of argument might explain the significant anisotropic isolation-by-distance observed on the North/South axis in France.

First, the significant anisotropic isolation-by-distance could be explained by the wind dispersal. Dispersion phenomena caused by wind currents are already established for *Culicoides*, mainly over the seas (Alba *et al*, 2004; Ducheyne *et al*, 2007; Eagles *et al*, 2012; Gloster *et al*, 2008; Hendrickx *et al*, 2008; Jacquet *et al*, 2016a; Murray, 1987; Murray and Kirkland, 1995). The reflecting dispersal events can be described as passive and active because recapture beam has been carried out downwind and upwind of the prevailing wind direction (Kluiters *et al*, 2015b). The map of the average wind direction in France over the last 10 years shows potential answers. It can be seen that the southernmost sampling sites (with the most dissimilar individuals) are in an area where the wind direction is different from the rest of France. It should be noted that the diffusion of the BTV, and thus the dispersion of the *Culicoides*, has already been associated with the wind direction. For example, 2% of BTV infections occurred at distances greater than 31 km (Sedda *et al*, 2012; Sedda *et al*, 2015) during the 2006 epizootic. The South of France was initially sampled, however a poor conservation of *Culicoides* DNA did not allow us to perform barcoding and microsatellite genotyping of these populations. The study of this geographical area could thus complete this study and potentially identify a stronger genetic structure. This will make it possible to decide whether this premise of differentiation observed in this study is due to an older phylogeographic structure, a different wind dispersion in this geographical area, or a random pattern. It is also questionable whether the very low indications of genetic structure could be due to differences in genetic diversity between the south and the north. A lower genetic diversity in the north could then be a sign of an expansion process that would contrast with the interpretation of a large-scale dispersion (Dlugosch and Parker, 2008; White *et al*, 2014). Indeed the spatial expansion of populations is generally accompanied by a gradients of reduction in the genetic diversity of the population during the expansion process, caused by serial founder events, creating a genetic bottleneck through a founding effect (Bock *et al*, 2015; Girod *et al*, 2011). Population sampling at the same scale could allow comparison and estimation of genetic diversity.

Second, anisotropic isolation-by-distance may be due to an artifact due to sampling methods. Indeed, if the extent of the sampling varies depending on the directions, the distances projected from the angles may represent different distance distributions and lead to the over-
representation of strong genetic metrics value in the direction of the scatterplot and the hit of
positive correlation signals. Moreover, the absence of correlation does not necessarily mean
more gene flow but an absence of isolation-by-distance which can also result from a strong drift
and thus less gene flow, a drift which depends on both dispersion and population sizes. It is
therefore essential to use as a complement, as we have done here and come up with a similar
result, an approach that weights the geographical distances between populations according to
their orientation with respect to a given angle axis passing through the barycenter, as
implemented in PASSAGE.

Consideration for future work

Our results underline the importance of methodological development of wind dispersal models
of *Culicoides* on land and not only over water. However, studies of landscape genetics remain
indispensable and complementary in order to improve the accuracy of predictive models for
*Culicoides* dispersal over land through integration of meteorological, landscape and activity-
based parameters previously tested and validated in the landscape genetics work flow. In
addition, population sampling of the same species at European level, *i.e.* a very large proportion
of the known range of the species, would be necessary to observe a more marked structuring at
European level and to estimate more precisely the gene flow. For example, *C. imicola* in Europe
is structured in two large genetic clusters, “Europe central cluster” and “Europe western cluster”
(Jacquet *et al*, 2016b). Moreover, the question of the phylogeographic history of *C. obsoletus*
is still very little explored in view of its geographical distribution throughout Europe and North
America. On the contrary, the study of the phylogeography of *C. imicola*, the main Afrotropical
vector, has shown its distribution range out from the northern part of sub-Saharan Africa to the
Mediterranean basin (Jacquet *et al*, 2015). This type of study would make it possible to estimate
the effective size of the populations, a key factor in the dispersal of *C. obsoletus*. In future study,
the use of High-throughput sequencing approaches using markers such as ddRadseq (Double-
digest restriction-site associated DNA sequencing) can provide greater resolution in view of the
large number of SNPs revealed (single nucleotide polymorphism) at a local scale and improve
our understanding of the active and passive dispersal of *Culicoides*. It could also be relevant to
look at on scale with finer genetic information. This could include more microsatellite markers
or SNP to improve genetic resolution and observe the matching and assignment of each
individual.
**Conclusions**

For the first time, this study provided the first complete landscape genetic analysis of *C. obsoletus*, a major vector species of animal viruses in Europe. This study showed that the large scale dispersion is able to homogenize the genetic structure of populations at the scale of a country. Our results demonstrate, for *C. obsoletus*, a very important inland dispersal and vectorization capacity, which has to be considered for further works on vector competence and epidemiological modeling of disease transmission. Wind direction could be a key factor in the dispersal of many vector insect species. Futures studies should increase their geographical extent to cover the entire species distribution area and better apprehend the limits of *Culicoides* gene flow. In addition to this biological knowledge, the study highlights several important areas for methodological improvement that may currently limit the inclusion of wind direction in landscape genetic analyses.

**Abbreviations**

BTV: bluetongue virus; SBV: Schmallenberg virus; *cox1*: cytochrome c oxidase subunit 1; N: number of individuals. SSR: simple sequence repeat. SNP: single nucleotide polymorphisms.

**Acknowledgement**

SD and MG are supported by the *Fonds National de la Recherche Scientifique* (FNRS, Belgium). AM and MJ have received a fellowships from the H2020-727393 PALE-Blu project funded the EU. We are most grateful to Hélène Vignes, Ronan Rivallan and Xavier Argout at the GPTR (Grand plateau technique regional de génotypage) platform in Montpellier, France for their help in development of microsatellite markers and Aurore Manez for her help in genotyping. We would like to thank Sylvain Piry, Karine Berthier and Marie-Pierre Chapuis, for their constructive exchanges which helped to improve the relevance of our conclusions. The authors are grateful to all partners who assisted in *Culicoides* sampling and shipping. We thank William Wint (University of Oxford, PALE-Blu partner) for his continued support in extracting and preparing meteorological data.
**Ethics approval and consent to participate**

Not applicable.

**Consent for publication**

Not applicable.

**Availability of data and materials**

All data generated or analyzed during this study are included in this published article and its additional files. The newly generated sequences were submitted in the GenBank database under the accession numbers MT828832-MT828844. All cox1 sequences are available on request.

**Competing interests**

The authors declare that they have no competing interests.

**Funding**

This study was partially funded by the EU grant H2020-727393 PALE-Blu and by the VectorNet project (OC/EFSA/AHAW/2013/02-FWC1) funded by the European Centre for Disease Prevention and Control (ECDC) and the European Food Safety Authority (EFSA). The contents of this publication are the sole responsibility of the authors and do not necessarily reflect the views of the European Commission, of the ECDC or of the EFSA.

**Author contributions**

AM, KH and CG designed the study. TB contributed to gather *Culicoides* samples and IR, MD contributed the identification of *Culicoides*. AM, LG, and M de W have performed molecular biology manipulations. AM, SD and MJ, developed techniques and workflow for landscape genetics analysis. AM analyzed the data. KH and CG contributed to the manuscript firstly written by AM. All authors read and commented the final manuscript version.