Introduction

1Ce jeu de données, publié dans trois articles (De Meeûs et al., 2002a, 2004a, 2004b), représente un excellent exercice, car nous allons y rencontrer bon nombre de situations décrites dans le chapitre précédent. Nous allons entièrement décortiquer une nouvelle fois ce jeu de données avec les mêmes méthodes, mais aussi avec des outils plus récents que ceux qui avaient été utilisés à l’époque, ce qui sera aussi intéressant. Nous repartirons de zéro en feignant d’ignorer ce qui a déjà été fait, comme s’il s’agissait d’un jeu de données non analysé. Le jeu de données complet est téléchargeable sur mon site web.

État des lieux

2Les tiques sont des acariens hématophages qui, au cours de leur repas sanguin, peuvent transmettre des maladies à leurs hôtes vertébrés. Dans l’hémisphère nord, ce sont elles qui sont responsables de la très grande majorité des maladies à vecteur des humains et, en particulier, de la transmission de la maladie de Lyme dont l’impact économique et en santé publique est reconnu (Gubler, 1998). Encore aujourd’hui, beaucoup reste à faire pour mieux comprendre l’épidémiologie de cette maladie et la variabilité des manifestations cliniques qui la caractérise (Hubbard et al., 1998). Les tiques sont typiquement des organismes difficiles à suivre sur le terrain, et des approches par marqueur moléculaire semblent donc pertinentes dans ce cas de figure. Après une tentative peu fructueuse avec les allozymes, avec seulement deux loci peu polymorphes (Delaye et al., 1997), des microsatellites ont été développés (Delaye et al., 1998). Seuls cinq loci polymorphes avaient pu être mis au point à l’époque, ce qui était vraiment peu. Nous allons ensemble voir que, malgré cela et les problèmes rencontrés, on peut quand même recueillir beaucoup d’informations pertinentes à l’aide des méthodes décrites dans ce manuel.

3En téléchargeant le fichier “IRTotBrut.txt”, vous aurez les données brutes obtenues sur des tiques adultes échantillonnées sur la végétation (donc non gorgées), sauf pour la Tunisie où les tiques étaient fixées sur des vaches. Le fichier “IRTotBrut.txt” est un fichier texte mais que l’on peut ouvrir sous Excel si on le souhaite. Le tableau 6 donne un extrait du fichier de données brutes. Le fichier comprend neuf colonnes. La première colonne donne le nom des sites où les tiques ont été échantillonnées. Il y a huit sites en Suisse (fig. 15) et un site en Tunisie. La deuxième colonne correspond à l’année d’échantillonnage, car certains sites ont été prélevés aux printemps 1995 et 1996 et d’autres uniquement au printemps 1996. La troisième colonne correspond au sexe de la tique (F pour femelle et M pour mâle). La quatrième colonne donne le nom codé des différents individus tiques. Ce codage individuel peut être utile si on fait des analyses individus centrées telle qu’une AFC ou une construction d’arbre sur distances interindividuelles. Enfin, les cinq dernières colonnes correspondent aux génotypes (en taille d’allèles) aux cinq loci microsatellites polymorphes définis dans Delaye et al. (1998).

4Il est à noter que les mâles sont tous homozygotes pour le locus IR08 car celui-ci s’est avéré lié à l’X, ce que nous ignorions lors de ce travail en 1996. Ils ont néanmoins été codés haploïdes chez les mâles dans le jeu de données téléchargeable. N’oubliez donc pas de les recoder homozygotes si vous souhaitez refaire les analyses dans les conditions initiales. Mais comme vous le constaterez, je préfère maintenant ne pas rendre homozygote ce qui ne l’est pas.

Figure 15

Premier recodage des données

5Certains logiciels n’aiment pas les noms longs et encore moins les accents ou autres signes cabalistiques. Par ailleurs, il est plus commode pour la lisibilité que tous les noms d’un même niveau aient le même nombre de caractères (alignement des colonnes). C’est pourquoi j’ai choisi de recoder dans IRTotBrut1.txt le nom des sites qui a été raccourci. Dans les données initiales, certains individus sont apparus avec trois ou quatre bandes à certains loci. Nous avions codé ces génotypes 333000 et 444000 pour les génotypes à trois et quatre bandes respectivement. Il convient de recoder ces données en données manquantes (000000). Nous reviendrons sur ces génotypes bizarres un peu plus tard, car ils s’avéreront utiles pour discuter des résultats des analyses de pedigrees. Nous allons procéder à une première analyse avec tous les échantillons afin de tester la panmixie locale et les déséquilibres de liaison entre loci. Nous allons pour ce faire créer un nouveau fichier où les sites et les dates seront distingués, mais aussi le sexe des tiques car on ne sait jamais à l’avance si des différences peuvent exister entre les deux sexes (Prugnolle et De Meeûs, 2002 ; Prugnolle et al., 2003), auquel cas les résultats obtenus pourraient s’en ressentir, mais surtout la discussion serait réorientée. Donc autant distinguer le sexe des individus dès le départ, quitte à ignorer ce facteur par la suite si on ne voit rien. Nous allons nommer ce fichier “IRTotTestPanmix.dat” et le mettre au format Fstat qu’il faut donc télécharger et ouvrir pour voir comment constituer un fichier à ce format. Vous pourrez aussi créer un fichier contenant le nom des sous-échantillons “IRTotTestPanmix.lab”, car un fichier de données Fstat ne contient que des chiffres. Ce fichier est constitué d’une colonne avec le nom des sous-échantillons. Vous pourrez aussi coder les données au format CREATE (Coombs et al., 2008) (qui n’existait pas au moment de réanalyser ces données en 2007) et vous servir de ce logiciel pour convertir ce fichier au format approprié.

Premières analyses : indépendance entre allèles dans et entre loci dans les sous-échantillons

6Nous allons donc tester s’il existe des déficits en hétérozygotes et des déséquilibres de liaison. Pour ce faire, il faut ouvrir Fstat. Une fois dans Fstat, il faut ouvrir le fichier “IRTotTestPanmix.dat” et cocher les cases qui vont nous être utiles ici (fig. 16). Si vous souhaitez voir apparaître les noms des sous-échantillons, il faut le spécifier par le menu “Options” de Fstat (cf. le premier recodage des données du chapitre 2 de cette deuxième partie pour une prise en main pas à pas de Create).

Figure 16

7Nous n’effectuons pas d’autres analyses pour le moment, car ces dernières pourraient être remises en cause par les résultats obtenus ici.

8La procédure de test de déséquilibre de liaison est assez lente, donc, si vous souhaitez que votre analyse finisse avant l’âge de la retraite, il vaut mieux dans tous les cas s’en tenir à l’option 5/100 pour le “Nominal level for multiple testing”. Mon ordinateur portable, dont l’horloge à 2.13 GHz et la mémoire vive à 2 Go témoignent d’une performance somme toute raisonnable, a mis quand même quatre heures pour effectuer cette première analyse dont le résultat est consultable dans le fichier “IRTotTestPanmix.out”. Si vous téléchargez Fstat 2.9.4 dans mon site web, vous pourrez choisir le nombre de permutations (je préconise 10 000). Que pouvons-nous voir dans ce fichier ?

9Les premières lignes donnent les fréquences des allèles pour chaque locus et chaque sous-échantillon, ainsi que sur l’ensemble (moyennes pondérée, W, et non pondérée, UW). Nous pouvons constater à cette occasion que chaque locus, sauf IR27, possède un très grand nombre d’allèles dont la plupart ne suivent en rien le modèle de mutation attendu de deux pas par deux pas (ce sont tous des dinucléotides). Dans ce cas, la plupart des allèles proviennent de mutations intervenues en dehors du motif microsatellite, dans les séquences flanquantes. Ce n’est pas dramatique même si non idéal, car cela favorise la survenue de problèmes liés au stuttering. Suivent les estimateurs de Nei, en particulier ceux des diversités géniques intra-sous-échantillons (H_s) et globale (H_T). Ensuite, les résultats des tests de déséquilibre de liaison sont donnés par paire de loci et par sous-échantillon et sur l’ensemble des sous-échantillons (mais toujours par paire de loci). La mention “Adjusted P-value for 5 % nominal level is : 0,000208” ne doit pas vous inquiéter. C’est le calcul du seuil de Bonferroni sur l’ensemble des tests réalisés. Comme il y a 24 sous-échantillons, cinq loci et donc 5(5 – 1)/2 paires de loci, cela correspond à 240 tests. Le seuil corrigé par la procédure de Bonferroni à α = 0,05 est donc α’ = 0,05/240 = 0,000208, seuil rarement (jamais ?) accessible, ce qui illustre une discussion que nous avons déjà eue précédemment. De toutes manières, nous ne regarderons ici que les tests multi-sous-échantillons (colonne “All”) et donc au pire, le seuil est à diviser par 10, ce qui est inutile puisque nous pouvons aussi constater qu’aucun déséquilibre de liaison n’est significatif. Les loci sont donc raisonnablement indépendants statistiquement les uns des autres. Nous pouvons donc sereinement oublier ces derniers et passer à la suite.

10Suivent les estimateurs de Weir et Cockerham dont un seul nous intéresse pour le moment, f, l’estimateur du F_IS, par locus, par allèle et sur l’ensemble des allèles, sur l’ensemble des loci. Puis suivent les résultats des jackknives et bootstraps et enfin des permutations. En compilant ces résultats dans le tableau 7 et la figure 17, nous constatons de très forts et très variables déficits en hétérozygotes (tous très significatifs avec des P-values toutes inférieures à 0,0001, visibles en fin de fichier).

Figure 17

11Ces fortes valeurs sont aberrantes étant donné qu’on sait qu’I. ricinus pratique une reproduction bi-parentale obligatoire. Des croisements systématiques entre apparentés pourraient-ils expliquer un F_IS = 0,5 ? Dans la réponse 11, on décrit comment obtenir une estimation grossière du taux de croisements frère-sœur b nécessaires pour expliquer un F_IS donné :

12Par conséquent, nous avons besoin ici de 4/5, soit 80 % de croisements frère-sœur pour expliquer nos données, ce qui est possible mais semble peu réaliste. Ixodes ricinus est en effet une tique triphasique qui change d’hôte pour chaque stade. Les adultes dont nous analysons la variabilité génétique ont donc subi deux phases de dispersion par des hôtes différents. Pour permettre un taux de 80 % de croisements frère-sœur, il faut admettre que 80 % des individus d’une même ponte restent ensemble au cours des différents stades (larvaire, nymphal et adulte) de leur vie.

13Il se pourrait, contrairement à ce qui est observé en laboratoire où aucun œuf non fécondé n’a pu éclore, que cette espèce pratique une parthénogénèse automictique d’un type qui augmente l’homozygotie (pour des descriptions des différents modes d’automixie, voir par exemple De Meeûs et al., 2007b). Seules les femelles sont en général capables de parthénogénèse. Il existe cependant une espèce de cyprès et une espèce de phasme où les mâles se reproduisent asexuellement (voir encore De Meeûs et al., 2007b) et une espèce de fourmi où mâles et femelles sont clonaux chacun de leur côté (Foucaud et al., 2010). Mais ce sont des exceptions. Si parthénogenèse il y a, les femelles devraient donc présenter de beaucoup plus gros déficits en hétérozygotes que les mâles (tous issus d’une reproduction croisée). Nous allons donc réanalyser le fichier en demandant à Fstat de nous donner les F_IS par sous-échantillon, puisque nous avons fort judicieusement, il faut bien l’avouer maintenant, d’entrée de jeu distingué les deux sexes.

14Sous Fstat vous ouvrez le même fichier “IRTotTestPanmix.dat” et vous décochez toutes les cases et cochez celle qui indique “Fis” dans le cadre “Per locus and sample statistics” comme indiqué dans la figure 18. Si vous souhaitez repérer encore une fois les noms des sous-échantillons, n’oubliez pas de signaler à nouveau l’existence du fichier “IRTotTestPanmix.lab” dans le menu “Options”.

Figure 18

15Quand vous lancerez “Run”, Fstat ouvrira une boîte de dialogue avec laquelle vous pouvez décider d’écrire les résultats de cette analyse dans un nouveau fichier. Dans le cas contraire, et c’est le choix que j’ai fait, le programme écrira les résultats dans “IRTotTestPanmix.out” à la suite des analyses précédentes (fin du fichier). Qu’y découvrons-nous ? Tout d’abord que Fstat tronque les labels plus longs que six caractères. Ce n’est pas grave, car nous avons toujours le même ordre Femelles Mâles pour chaque échantillon. Et puis il suffit (sous Excel c’est facile) de faire un copier-collage spécial/transposition à partir du fichier “IRTotTestPanmix.lab”. Ensuite, comme représenté dans la figure 19, construite à partir du fichier de sortie, nous pouvons voir, qu’à part pour le locus IR08, aucune tendance claire n’apparaît. Tous ces loci présentent des déficits importants et relativement variables, mais sans lien réel avec le sexe des tiques. Ce seraient plutôt les mâles qui auraient une tendance à présenter des déficits d’hétérozygotes plus importants (nous verrons plus loin une explication possible). Pour le locus IR08 par contre, avec un F_IS = 1 pour les mâles, il apparaît clairement que ce locus est situé sur le chromosome X et qu’il est donc haploïde chez les individus mâles.

Figure 19

16En fait, pour être précis, le locus IR08 avait été trouvé hétérozygote pour quatre individus mâles sur l’ensemble du jeu de données. Même si cela pouvait refléter des duplications toujours possibles (comme évoqué p. 124), nous avons choisi d’éliminer ces individus, car ils pouvaient correspondre à des erreurs de manipulations.

17Quoi qu’il en soit, il va donc falloir recoder les données à ce locus. Pour l’analyse des F_IS, les mâles devront en effet être codés en données manquantes (000000) au locus IR08. Nous allons donc créer un nouveau fichier “IRTotTestPanmixMalManqIR08.dat” à partir du précédent et refaire l’analyse globale du F_IS. Celle des déséquilibres de liaison, qui est un test génotypique, n’a aucune raison d’avoir été affectée par ce phénomène. Dans Fstat, nous cocherons donc les mêmes cases qu’en figure 16, à l’exception de celles concernant les déséquilibres de liaison.

18Dans le fichier de sortie “IRTotTestPanmixMalManqIR08.out”, nous constatons l’image suivante (voir aussi la figure 20) : rien ne change sauf pour le locus IR08 qui montre les plus basses valeurs de F_IS, mais qui restent très significativement (toutes les P-values sont inférieures ou égales au minimum possible 0,0001) au-dessus de la valeur nulle attendue sous panmixie. Notez au passage que je ne me sers des intervalles de confiance que pour illustration. Le F_IS global reste donc très élevé (0,39), inexplicablement variable entre loci et fort variable d’un site à l’autre. Ceci suggère un rôle possible pour des allèles nuls ou de dominance d’allèles courts. L’étape suivante sera donc de mettre en évidence l’existence de tels allèles et/ou de phénomène de dominance des allèles les plus courts.

Figure 20

Recherche d’allèles nuls et de dominance d’allèles courts

19Nous allons pour ce faire utiliser deux nouveaux logiciels. Micro-Checker va nous permettre d’estimer la fréquence des allèles nuls susceptibles d’expliquer, dans chaque sous-échantillon et pour chaque locus, les déficits en hétérozygotes observés. Micro-Checker permet également d’estimer si les données sont compatibles avec un bégaiement de la polymérase (stuttering) et/ou une dominance des allèles les plus courts. Pour la dominance des allèles courts, nous utiliserons également une méthode plus puissante que celle implémentée par Micro-Checker. Nous allons procéder à une régression généralisée pour la mise en œuvre de laquelle nous utiliserons le logiciel R (R-Core-Team, 2020 : voir la référence complète dans la bibliographie).

Convertir le fichier pour Micro-Checker et ouverture du logiciel

20Pour commencer avec Micro-Checker, nous avons besoin de transformer nos données au format Genepop qui est compatible avec ce logiciel. Ensuite, nous allons devoir créer un fichier spécial pour les données du locus IR08, lié au sexe, sans les mâles car sinon Micro-Checker risque de goûter moyennement la saveur de cette plaisanterie. Créons donc un fichier “IR08AllFem.txt” avec les données femelles pour le seul locus IR08 et un fichier “IRAutosomAll.txt” pour le reste des données. Attention, le fichier doit suivre des règles strictes sinon Micro-Checker refusera d’analyser les données. Référez-vous au fichier exemple fourni avec le logiciel et respectez les espaces et tabulations de la façon la plus scrupuleuse (ou utilisez Create).

21Lancez Micro-Checker et ouvrez “IRAutosomAll.txt” avec le menu “File”. Si tout se passe bien, vous observez l’ouverture de votre fichier avec vos données et différents menus et boutons en bas de l’écran.

Analyses des loci autosomiques du premier sous-échantillon par Micro-Checker

22Il y a un encadré en bas à gauche où il faut choisir le motif de chaque locus microsatellite. Il affiche par défaut le premier des loci (ici IR25) et un blanc pour le motif. Choisissez le motif “Mononucleotide” comme sur la figure 21.

Figure 21

23Nous avons déjà remarqué que nos loci microsatellites étaient peu orthodoxes. L’option mononucléotidique correspond en fait à l’option qui permet de faire face à toutes les situations. Cliquez ensuite sur le bouton “All” pour signaler que cette option est valable pour tous les loci, puis aller au locus IR27 pour le remettre en dinucléotide. Cliquez ensuite sur le bouton “Analyse” (un peu plus à droite). Apparaît alors une fenêtre d’avertissement comme celle présentée en figure 22. Comme il y a des données manquantes, Micro-Checker vous demande s’il faut ou non en tenir compte. Autrement dit, les données manquantes correspondent-elles à des homozygotes nuls (blancs) et faut-il les utiliser pour le calcul des fréquences des allèles nuls par la seconde méthode de Brookfield (1996) ? La réponse étant positive, cliquez donc directement sur “Proceed” sans vous poser plus de questions.

Figure 22

24Micro-Checker effectue plusieurs calculs et vous présente des résultats concernant le premier locus. Allez dans le menu “Tools” à “Nulls across loci” comme dans la figure 23 pour obtenir le tableau des fréquences de nuls dans le premier sous-échantillon, estimées selon différentes méthodes. Sélectionnez ce tableau avec la souris, copiez-le et sauvez-le dans un fichier (Excel, par exemple). Ensuite, regardez dans l’encadré en bas à droite (fig. 23) si le locus correspondant montre un problème de stuttering ou une dominance d’allèle court (« large allele dropout »). Si oui, notez-le dans le tableau que vous venez de créer pour sauvegarder les résultats de cette analyse puis, par le menu “Window’” (fig. 23) sélectionnez le locus suivant, etc. Vous constaterez qu’aucun locus ne présente de « stuttering » ni de dominance d’allèle court dans ce premier sous-échantillon.

Analyses des autres sous-échantillons, des autres loci autosomiques et du locus IR08

25Au centre et en bas, cliquez sur le bouton “Next Population” (voir fig. 23) pour analyser le sous-échantillon suivant en reprenant les mêmes étapes décrites en p. 131-133, jusqu’au dernier sous-échantillon. N’oubliez pas de copier le tableau des fréquences d’allèles nuls à chaque fois (dans le menu “Tools” à “Nulls across loci”, fig. 23). Ensuite, vous ferez la même chose pour le locus lié au sexe, IR08, en ouvrant le fichier correspondant “IR08AllFem.txt”.

Figure 23

Bilan des analyses avec Micro-Checker

26Nous avons constitué un fichier de résultats avec les fréquences d’allèles nuls probables, l’existence ou non de stuttering et de dominance d’allèles courts. Nous ne gardons que la méthode 2 de Brookfield (1996) qui tient compte des données manquantes (blancs) comme des homozygotes nul/nul. Dans ce fichier, nous allons également insérer le nombre d’individus génotypés pour chaque locus (copiés à partir des fichiers de sortie Fstat), la fréquence attendue sous panmixie (fréquence précédente au carré) des allèles nuls pour chaque locus dans chaque sous-échantillon et sur l’ensemble des sous-échantillons, le nombre de blancs observés (compter les 000000 dans chaque sous-échantillon et sur l’ensemble), l’effectif corrigé (individus génotypés + blancs) et enfin le nombre de blancs attendus sous la double hypothèse qu’il y a panmixie et que les allèles nuls expliquent les F_IS en totalité. Le tableau 8 donne un aperçu du fichier final pour le locus IR08.

27Pour vérifier que ces résultats expliquent correctement les F_IS observés, on peut comparer la proportion de blancs observés avec celle attendue sous l’hypothèse que les allèles nuls expliquent la totalité de ces F_IS. Un test binomial unilatéral avec comme fréquence attendue p_B2², un nombre de réussite égal aux blancs observés pour un nombre d’essais de N’, semble ici approprié. On préfère ici un test unilatéral, car ce qui nous intéresse est de savoir si on a oui ou non moins de blancs qu’attendus. On peut facilement effectuer ce test sous R.

28Il nous faut donc lancer R et dans la fenêtre de commande taper l’instruction :

29binom.test(Blancs observés, N’, p = p_B2², alternative = “less”)

30Pour des raisons de recherche de puissance et pour limiter le nombre de tests dont la multiplication est toujours problématique (voir p. 84 en première partie), on ne fera les tests qu’avec les valeurs totales pour chaque locus. Pour le locus IR08, cela correspond aux valeurs de la dernière ligne du tableau 8. Pour ce locus, la commande devient donc :

31binom.test(4, 186, 0.0253, alternative="less")

32Faites bien attention de respecter strictement le format (en particulier, les majuscules et minuscules sont reconnues comme des caractères différents sous R). Ici “less” signifie que le test est unilatéral dans le sens des plus petites valeurs (H1 : il y a moins de blancs observés qu’attendus) (l’instruction devient “two.sided” pour un bilatéral et “greater” pour l’autre test unilatéral). Une fois que vous avez tapé cette instruction dans R, tapez sur la touche “Entrée” et le test se fait. La P-value du test est, pour IR08, non significative (P-value = 0,4919). Les allèles nuls sont donc bien suffisants pour expliquer les déficits en hétérozygotes observés à ce locus chez les femelles, d’autant plus qu’il semble aussi exister des phénomènes de stuttering à ce locus. Pour les autres loci, on procède de la même façon. On trouve ainsi que pour les loci IR 25, IR27 et IR32, la fréquence des blancs observés est significativement inférieure à celle des blancs attendus si les allèles nuls devaient expliquer les déficits en hétérozygotes. C’est un problème car, par un phénomène de cercle vicieux, moins les allèles nuls expliquent un déficit en hétérozygotes, moins le nombre de blancs observés correspond aux attendus. Pourquoi cela ? Simplement parce que si on attend naturellement plus d’homozygotes en général, alors on devrait observer encore plus d’homozygotes nuls (blancs), en particulier (ce raisonnement ne marche cependant pas très bien s’il s’agit d’un effet Wahlund). Par ailleurs, la variance entre loci ainsi que le fait que les nuls expliquent très bien les déficits observés pour IR08 (voir plus haut), mais aussi pour IR39 (P-value = 0,312) pourraient nous inciter à exclure des causes biologiques du type régime de reproduction ou effet wahlund (voir plus loin). Notons que des phénomènes de stuttering ont été détectés pour IR25, mais seulement dans deux sous-échantillons. Pour IR32 et IR27, Micro-Checker n’a pas détecté ce phénomène pas plus qu’il n’a détecté de dominance d’allèles courts. Cependant, Micro-Checker ne travaille que dans chaque sous-échantillon de façon isolée, ce qui peut représenter une forte perte de puissance. Dans le paragraphe qui suit, nous allons utiliser une autre technique pour détecter d’éventuelles dominances d’allèles courts.

Détection de dominance d’allèles courts par la méthode de régression multiple

33Pour ce faire, nous aurons besoin de connaître, pour chaque locus et dans chaque sous-échantillon, la valeur du F_IS pour chaque allèle. On peut demander à Genetix de le faire en choisissant à chaque traitement le locus et le sous-échantillon à analyser, en n’oubliant pas de zapper les mâles au locus IR08. On peut aussi créer autant de fichiers Fstat qu’il y a de sous-échantillons à analyser, ensuite, et parce que malheureusement Fstat ne permet pas d’analyser qu’un seul sous-échantillon, il faut créer dans chaque fichier une deuxième population fictive, de taille identique à celle à analyser et fixée à tous les loci (par exemple, tous homozygotes 170170, 150150, 123123, 235235, 129129 pour les cinq loci respectivement). Il s’agit ensuite de récupérer dans chaque sous-population les F_IS de chaque allèle pour chacun des cinq loci et de créer cinq fichiers de données (un par locus) contenant pour chaque allèle son F_IS, sa taille (on s’en doute), le sous-échantillon, sa fréquence allélique p dans ce sous-échantillon, le produit p(1 – p), le nombre d’individus génotypés dans ce sous-échantillon N et enfin le produit p(1 – p)N. Le tableau 9 donne une idée de la forme de ce fichier pour le locus IR08 que j’ai appelé “IRTotL08MalManqFisAllSizeL08.txt”. Pour fabriquer ce fichier, une feuille de calcul Excel est idéale, ensuite il suffit d’enregistrer le fichier en format texte seul.

34On peut aussi utiliser Genetix qui permet l’analyse d’un seul sous-échantillon, mais dont les sorties sont moins commodes à importer dans Excel (à vous de voir).

35La colonne Npq, qui donne en fait le résultat du produit Np(1 – p), nous servira à pondérer notre régression par la taille des échantillons, mais en donnant aussi plus de poids aux allèles de fréquences proches de 0,5 (les plus polymorphes). On fait les mêmes fichiers avec les quatre autres loci. Nous allons maintenant analyser ces données avec le logiciel R.

36Ouvrez R et dans le menu “Fichier” cliquez dans “Changer le répertoire courant...”, et allez dans le répertoire où vous avez stocké vos fichiers de données. Dans la console de travail de R, tapez la suite de commandes, chacune suivie d’un retour chariot (touche “Entrée”) :

37> data<-read.table("IRTotL08MalManqFisAllSizeL08.txt", header=TRUE)

38qui signifie que le tableau de données “data” est contenu dans le fichier nommé et que la première ligne contient le nom des colonnes. N’oubliez pas que les données manquantes se notent “NA” en majuscules et non “000000”.

39> attach(data)

40qui signifie que ce tableau doit être chargé en mémoire1.

41> loc8<-glm(data, formula = Fis ~ poly(Allele, 2) + Sample + Year, family = gaussian, weights = Npq)

42où loc8 est le nom d’un modèle linéaire généralisé utilisant le tableau “data” et dont la régression tente d’expliquer la valeur du F_IS en fonction de la taille des allèles selon un polynôme d’ordre 2 ou quadratique (qui s’est avérée plus proche de ce qui se passe dans le cas qui nous intéresse), du sous-échantillon d’origine et de l’année. Le sexe n’a ici aucune importance puisqu’il n’y a que des femelles. Nous ne testons l’effet d’aucune interaction entre variable, car en fait je ne vois aucune raison pour qu’il en existe. Pensez à respecter les majuscules s’il y en a, car R les reconnaît comme telles. Tapez enfin :

43> anova(loc8, test="F")

44qui renvoie à une analyse de variance utilisant la statistique F (se référer à un livre de statistique pour approfondir ces notions) et donne le résultat suivant :

45Ici, on voit que seul le sous-échantillon influence la valeur du F_IS (allèles nuls, stuttering variable dans l’espace ?) qui n’explique que 14,63 % de la dispersion (100 × 3,1604/21,616), tout en étant très significatif.

46On utilise un test F, car on a supposé que la distribution des F_IS suit plus ou moins une courbe de Gauss (données continues en cloche symétrique), ce qui est sûrement inexact mais ne risque guère de modifier le résultat dans un sens dramatique.

47Pour les loci suivants, nous aurons besoin de distinguer le sexe des tiques.

48Avec le locus IR25, l’analyse du fichier “IRTotSexSepFisAllSizeL25.txt” est la suivante :

49Ce qui aboutit au tableau de résultat :

50On aboutit à une conclusion similaire à la précédente, puisque ni le sexe ou l’année ni la taille des allèles ne comptent avec seulement environ 6,13 % de la déviance expliquée par le site qui est moins spectaculairement significatif que précédemment.

51Pour le locus IR27, le tableau obtenu est différent :

52En effet, comme nous pouvons le déduire du tableau ci-dessus, le site (Sample) explique 19,96 % de la dispersion des points (5.181/25.9549) et la taille des allèles (poly(Allele, 2)) en explique 15,86 % (4.1186/25.9549) et sont tous les deux très significatifs (souligné par les trois étoiles). Ils expliquent ainsi 35,83 % de la variance. Cette valeur est conséquente eu égard à l’importante variance résiduelle attendue en général pour un estimateur de statistique F. Comme le montre la courbe décrite dans la figure 24, la relation entre F_IS et taille des allèles est négative (si on exclut les trois premiers points, ce qui ne changerait rien eu égard aux intervalles de confiance), ce qui peut donc être interprété par une dominance des allèles les plus courts.

Figure 24

Les intervalles de confiance à 95 % ont été obtenus avec .
Pour ce faire, les singletons (tailles d’allèles présents une seule fois comme 112 et 131) ont été réunis à la classe la plus proche.

53Pour le locus IR32, on observe le résultat suivant :

54On voit qu’en plus du site, le sexe des tiques a un effet significatif, ce qui signifie que nous avons eu raison d’en tenir compte et nous verrons ensuite pourquoi.

55Pour le locus IR39, le tableau obtenu est le suivant :

56Le site joue une fois encore de façon significative, mais aussi la taille des allèles, même si cette dernière n’explique même pas 2 % de la déviance et est peu significative. Par ailleurs, la figure 25 montre que la relation (augmentation globale du F_IS avec la taille des allèles) n’est pas compatible avec une dominance des allèles courts. On peut donc attribuer ce résultat au hasard et au nombre de tests effectués qui augmente la probabilité d’obtenir quelque chose de significatif par hasard (revoir la première partie de ce manuel, p. 84). Rappelons que pour ce locus, les allèles nuls s’étaient avérés suffisants pour expliquer les déficits en hétérozygotes observés. Il est plus raisonnable ici de considérer ce résultat comme fortuit.

Figure 25

Bilan de l’analyse des déficits locaux en hétérozygotes

57Pour les loci IR08 et IR39, les allèles nuls semblent pouvoir expliquer les forts et variables F_IS observés. Pour IR27, les allèles nuls et la dominance des allèles courts offrent conjointement une explication satisfaisante. Seul le locus IR32 offre des déficits énormes et non expliqués par les allèles nuls, le « stuttering » ou la dominance des allèles courts. Cependant, sachant que le « stuttering » n’a pu être testé que sous-échantillon par sous-échantillon (manque de puissance), que la plupart des allèles se suivent à un pas sur ce locus et compte tenu de ce que nous trouvons aux autres loci, il est possible qu’ici aussi les déficits observés proviennent d’un problème technique.

58Je peux ajouter ici qu’un module (package) de R, appelé “R-Commander”, dont je n’ai appris l’existence qu’après la rédaction de ce chapitre, permet d’accéder aux analyses effectuées dans ce paragraphe à l’aide de menus déroulants plus conviviaux que le mode commande strict.

Recherche d’une structure cachée (effet Wahlund)

Introduction

59Nous allons dans un premier temps continuer de considérer les femelles et les mâles séparément. On sait en effet qu’il y a une structure génétique spécifique pour chaque sexe dans ce jeu de données. Même si nous analyserons ceci plus tard, il n’est pas inutile de poursuivre la recherche d’explications des déficits en hétérozygotes avant d’aborder cet aspect. Nous allons donc analyser tous les sous-échantillons (mâles et femelles séparées) pour obtenir l’information sur le plus grand nombre de réplicas possibles. Ensuite, nous nous concentrerons sur 1996 en réunissant les mâles et les femelles pour faire des tests.

60Le but du jeu sera ici d’utiliser l’information multilocus de chaque individu, dans chaque sous-échantillon afin de vérifier à l’aide du logiciel BAPS (voir p. 101-105 en première partie et le tableau A1 en annexe), si certains individus peuvent être regroupés sur la base de leur ressemblance génétique. BAPS va ensuite explorer de façon itérative et répétée, en suivant plusieurs chaînes de Markhov (ou une chaîne stochastique d’optimisation suivant les versions) afin de trouver la meilleure partition (celle qui regroupe le mieux les individus) dans le sous-échantillon analysé. La partition définit un nombre donné de clusters (sous-unités) composés chacun d’un certain nombre d’individus du sous-échantillon. La qualité d’une partition se définit par un savant calcul dans le détail duquel je serai bien incapable de rentrer, mais qui dépend de la distance génétique entre les groupes définis, par rapport aux autres partitions explorées durant le processus. Il est aussi expliqué, dans les articles décrivant le logiciel, qu’une hypothèse du modèle utilisé dans l’algorithme est que les « clusters » qui composent la partition sont en équilibre de Hardy-Weinberg. Je ne suis pas certain de bien comprendre ce qui est entendu par là dans la mesure où mon expérience m’a montré que la plupart des partitions obtenues ne sont pas conformes à cet équilibre, voire même en sont très éloignées. J’ai également pu observer cela avec STRUCTURE qui fait la même hypothèse. Comme discuté dans la première partie de ce manuel, beaucoup reste à explorer concernant le fonctionnement de ces méthodes dans différentes situations. Il faudra donc vérifier si la partition obtenue (car le logiciel en donne toujours une) correspond à quelque chose de viable et pas seulement une vue de l’esprit.

61Si la partition a réellement mis en évidence des groupes cryptiques au sein des sous-échantillons susceptibles d’expliquer en partie (effet Wahlund) nos fameux déficits en hétérozygotes, il faudra ensuite trouver et explorer les hypothèses susceptibles d’expliquer le plus raisonnablement possible (mais en aveugle) ces résultats (espèces ou races d’hôtes cryptiques, sous-structures familiales, isolement par la distance entre individus sur de courtes distances).

62Il existe d’autres logiciels qui en principe font la même chose. L’avantage de BAPS réside dans sa convivialité, dans le fait qu’il accepte des fichiers de type Genepop (un peu modifiés) et qu’il m’a toujours donné de bons résultats. Le logiciel STRUCTURE est par exemple beaucoup moins commode à utiliser (et c’est un euphémisme) et, sur un même jeu de données (glossines), n’a pas offert de partitions aussi satisfaisantes que BAPS (Ravel et al., 2007). Des études comparatives de différents logiciels de clustering sont en cours, mais la longueur et la quantité des analyses font que des résultats concrets ne seront sans doute pas disponibles avant la sortie du présent ouvrage. Vous verrez aussi l’application d’un autre logiciel de même nature, Flock, plus loin dans cette partie.

Construction des fichiers BAPS

63Il faut construire un fichier pour chaque sous-échantillon. Le type est semblable à un fichier Genepop, mais avec des tabulations comme dans la figure 26 (symbolisées par des →) qui donne un exemple pour le fichier des mâles de Staadswald. On note que les mâles sont codés homozygotes pour IR08 afin que l’information multilocus soit préservée pour cinq loci. Par commodité, j’ai appelé ce fichier “IRTotBrut1Stad95M.gen”, mais vous faites comme bon vous semble.

Figure 26

64Ensuite, il est commode de créer un fichier texte contenant le chiffre 30 répété un grand nombre de fois (ici 50 fois), avec un espace entre chaque répétition et sur une seule ligne. Le logiciel BAPS vous demandera en effet de taper un nombre maximal probable pour les clusters. Ici, 30 m’est apparu comme largement raisonnable compte tenu des tailles de sous-échantillons. C’est à partir de ce chiffre que BAPS démarre et recherche une partition la plus probable en se limitant à ce nombre maximum de clusters. Le logiciel reprend ensuite le processus autant de fois que l’on a rentré ce chiffre (ici 50) et ne gardera que la meilleure de toutes les partitions explorées. Avoir tapé 50 fois ce chiffre dans un fichier permet de copier et coller cette séquence directement sans avoir à la retaper pour toutes les analyses. J’ai appelé ce fichier “50fois30.txt” (quelle imagination !).

Analyse des fichiers par BAPS

65Vous avez bien entendu installé BAPS sur votre machine et créé tous les fichiers nécessaires (il y en 24 normalement). Il faut maintenant lancer BAPS en cliquant sur BAPS4_RUNME.EXE. Le logiciel ouvre deux fenêtres, une fenêtre Dos dont il n’est pas vraiment nécessaire de se préoccuper maintenant et une fenêtre d’interface type Windows avec des menus que nous allons utiliser. Il est important de commencer par créer un fichier résultat. Pour ce faire, cliquez sur “File”, “Output File” et “Set” et créez un fichier en tapant son nom et en le plaçant dans le répertoire qui vous convient le mieux (là où sont vos données) (voir fig. 27).

Figure 27

66Il vaut mieux garder un nom de fichier qui permette de retourner ensuite au fichier de données correspondantes. Ici, le premier fichier analysé sera “IRTotBrut1Ber95F.gen” (femelles de Berne 1995), je choisis donc ici de nommer et créer le fichier résultat “IRTotBrut1Ber95FBAPSRes.txt”. Ensuite, il faut cliquer sur le bouton “Clustering of individuals” (fig. 27). Apparaît alors une nouvelle fenêtre de dialogue qui vous propose différents formats de fichiers de données (fig. 28). Choisissez bien entendu le format Genepop en cliquant sur le bouton correspondant. Une fenêtre qui s’ouvre vous permet de naviguer vers le répertoire où se trouve IRTotBrut1Ber95F.gen que vous sélectionnez (soit en tapant son nom complet, soit en tapant *.gen et retour chariot et en double cliquant sur le fichier). Une question vous est alors posée “Do you wish to save pre-processed data?”, cliquez sur “No”.

Figure 28

67C’est alors qu’apparaît une petite fenêtre permettant de sélectionner le nombre maximum de clusters, ainsi que le nombre de chaînes d’itérations à effectuer (fig. 29), comme expliqué en p. 142. Supprimez le chiffre par défaut (20) et remplacez-le par la chaîne de 30 que vous copiez à partir de “50fois30.txt”, collez cette chaîne dans la case idoine et cliquez sur “OK”.

Figure 29

68Les calculs démarrent et se poursuivent jusqu’à la fin où la meilleure partition est sauvée dans “IRTotBrut1Ber95FBAPSRes.txt”. Apparaissent un graphique censé représenter la partition (clusters de différentes couleurs), dont on ne va pas se servir, ainsi qu’un dialogue final vous demandant si vous souhaitez sauver ces données en vue d’une analyse ultérieure. Répondez non. Ceux qui souhaitent plus de détails sur BAPS et ses différentes possibilités et menus sont invités à consulter la documentation livrée avec le logiciel.

69Il s’agit ensuite de répéter le processus avec chacun des sous-échantillons. Ensuite, on charge le jeu de données brutes afin de le modifier. N’oubliez pas de créer un nouveau fichier de sortie à chaque fois. Dans chaque fichier de résultat BAPS sont donnés les clusters avec les individus qu’ils contiennent. Ces individus sont identifiés par leur rang d’entrée dans le jeu de données (1,2,3...). Par exemple, pour les femelles de Berne 1995, le fichier de résultat donne (en début de fichier) :

70Le nom du fichier analysé est suivi de l’effectif de l’échantillon, du nombre de clusters dans la meilleure partition et de la valeur du Log de la valeur marginale de vraisemblance ou Log(MV) qui sert de critère à BAPS pour sélectionner la meilleure partition, c’est-à-dire celle qui présente le plus petit Log(MV). C’est bon à savoir si on souhaite relancer BAPS sur les mêmes données afin de voir s’il trouve une partition meilleure au deuxième essai. Enfin, la partition est donnée. Dans le jeu de données, il faut donc maintenant ajouter une colonne avec le numéro de cluster BAPS auquel chaque individu appartient. Il faut le faire pour tous les sous-échantillons (cf. tabl. 10). Attention, vous allez peut-être trouver des partitions légèrement différentes des miennes et avec des labels de clusters différents, c’est normal.

71Il faut ensuite créer un nouveau fichier de données où chaque sous-échantillon initial se retrouve subdivisé en autant de sous-échantillons que de clusters de BAPS qui le composent (12 pour les femelles de Berne 1995). Sous un éditeur quelconque vous fusionnez les colonnes 1, 2, 3 et 6 du tableau 10, ce qui donne pour la première ligne quelque chose du style Ber95F1. N’oubliez pas de trier les données pour que les clusters apparaissent dans l’ordre dans chaque sous-échantillon initial. Appelons le fichier contenant ces données modifiées “IRTotBAPSClustMalHomoMFSep.txt”.

72Ce n’est pas fini, car il faut maintenant coder en données manquantes le locus IR08 chez les tiques mâles. Rappelez-vous que, dans un souci de puissance, nous les avions artificiellement rendus homozygotes afin que les mâles soient pris en compte pour ce locus dans l’analyse BAPS. Maintenant, nous souhaitons calculer les nouveaux F_IS de cette partition afin de voir si elle chute par rapport au jeu de données initiales. Le génotype des mâles au locus IR08 doit donc en effet être recodé 000000, car ils ne doivent pas rentrer en ligne de compte dans le calcul du F_IS. Pour ce faire, il est commode soit de faire un petit programme (pour ceux qui savent), soit d’utiliser la fonction conditionnelle d’Excel. Il s’agit de créer une colonne sexe en A dans le jeu de données “IRTotBAPSClustMalHomoMFSep.txt” avec le sexe des individus (F ou M), dans une colonne libre (en H après IR39) on tape en ligne 2 (ligne du premier individu) :

73SI(A2="M";"000000";C2), ce qui aura pour effet d’écrire “000000” dans la case H2 si l’individu est mâle ou de recopier le génotype de la femelle au locus IR08 (contenu dans la case C2). On copie ensuite H2 et on le colle de H3 à H726 (normalement la fin du fichier). On sélectionne les cases H2 à H726, on les copie et on fait un collage spécial (on veut ne coller que la valeur et non la formule) sur C2. On supprime les colonnes H et A et on sauve en texte seul sous le nom “IRTotBAPSClustMalManqIR08MFSep.txt”. Supprimez aussi le label de la première colonne (c’est pour Genetix qui ne désire que le nom des loci).

74Nous allons maintenant recalculer les F_IS par locus et sur l’ensemble, avec intervalles de confiance. Nous sommes paresseux et pour ne pas avoir à supprimer les clusters d’un individu pour lequel le calcul ne se fera pas, et étant donné que Fstat ne prend pas plus que 200 sous-échantillons (avec mes partitions je me retrouve avec 368 sous-échantillons), nous allons importer notre nouveau fichier sous Genetix. Lancez Genetix et allez dans le menu “Fichier”, sélectionnez “Importer” et sélectionnez “Texte avec séparateur” et sélectionnez le fichier. Un menu apparaît et si vous avez fait comme moi, vous devez cocher les cases comme dans la figure 30. Quand cela est fait, cliquez “OK”. Si le fichier est correctement chargé, cliquez dans le menu “Fstats” et sélectionnez “Weir & Cockerham”. Cliquez OK dans la nouvelle fenêtre si vous ne changez pas le nom du fichier de sortie proposé “IRTotBAPSClustMalManqIR08MFSep.res”. Après il faut prendre les résultats dans ce fichier en prenant garde que ce qui est annoncé comme écart-type des jackknives correspond à l’erreur standard de Fstat. Il s’agit de comparer maintenant les F_IS de chaque loci et leurs intervalles de confiance de jackknife sur populations (voir p. 73-76 en partie 1) avant et après clusterisation par BAPS, ainsi que les valeurs globales et leur intervalle de confiance de bootstrap sur les loci (voir p. 73-76 en partie 1). La compilation des résultats prend alors la forme de ce qui est représenté dans la figure 31. Dans cette figure, il est aisé de voir que les clusters de BAPS présentent des déficits en hétérozygotes significativement inférieurs au F_IS de départ. Un test de rang de Wilcoxon pour données appariées confirme cela. Pour effectuer ce test sous R, il faut construire un fichier avec une colonne “Delta” où chaque ligne correspond à un locus.

Figure 30

75Chaque valeur représente la différence entre le F_IS brut et le F_IS BAPS au locus correspondant (ici cinq valeurs). Appelons ce fichier “DeltaFisBrutBAPS.txt”. Ensuite, sous R les commandes sont les suivantes :

76Le test est unilatéral, car ce que nous recherchons est bien un effet Wahlund. Nous attendons au départ une chute du F_IS, d’où l’instruction “greater”. La P-value = 0,031 obtenue est significative. Notons aussi que la plupart des loci, mis à part IR08, gardent un fort F_IS qui provient probablement des allèles nuls et autre dominance des allèles courts. Ces déficits restent très significativement au-dessus de 0 (fig. 31), ce qui rend bien compte du fait que “Hardy-Weinberg” n’est pas une nécessité pour parvenir à une partition. Par ailleurs, le F_IS fait mieux qu’être faible pour IR08, il est négatif, ce qui est effectivement ce que nous attendons chez une espèce dioïque pangamique.

Figure 31

77Il semble donc bien y avoir un effet Wahlund, contrairement à ce que la variance du F_IS entre loci pouvait laisser prévoir. Reste à déterminer si cet effet provient d’une micro-structuration (en groupes familiaux, par exemple) ou de la présence d’espèces (ou races d’hôtes, ou groupes adaptatifs ou écotypes) cryptiques. Afin d’essayer d’argumenter dans un sens ou l’autre, on peut essayer de regarder l’organisation de ces différents clusters. En principe, si on a à faire à différentes espèces, ces dernières devraient apparaître clairement. Si on effectue un arbre à partir d’une matrice de distance inter-clusters, ces derniers devraient être regroupés selon l’espèce à laquelle ils appartiennent en groupes séparés par des branches relativement longues comparées aux branches séparant chaque cluster (géographique, en principe) à l’intérieur de chaque espèce. Selon Takezaki et Nei (1996), la méthode du Neighbor-Joining (NJTree) sur distances de corde de Cavalli-Sforza et Edwards (1967) est une bonne solution. La matrice est obtenue en important “IRTotBAPSClustMalHomoMFSep.txt” dans Genetix2, en cliquant sur le menu “Distances” puis “Sur données réelles” et en sélectionnant “Cavalli-Sforza & Edwards”. On copie la matrice obtenue afin de l’incorporer dans un fichier de type MEGA (Kumar et al., 2004) pour matrice de distances (ouvrir le fichier “IRTotBAPSClustMalHomoForNJTREENmini3CSE.meg” avec un éditeur de texte pour voir un exemple). Afin de limiter le nombre de branches et le poids des clusters ne contenant qu’un seul ou deux individus, je n’ai gardé que les clusters d’au moins 3 individus. L’arbre obtenu n’en est pas plus lisible pour autant et ce qui en ressort, c’est que les plus longues branches sont toujours celles séparant les clusters sans que se dégage une quelconque hiérarchie (on parle de râteau). Ceci plaide davantage en faveur d’une micro-structuration locale forte avec une différenciation géographique faible. On peut alors recommencer l’ensemble des opérations (BAPS->Genetix->MEGA) sur les échantillons de 1996 seuls et en ne séparant pas les mâles des femelles. Sur l’arbre obtenu, on ne voit pas mieux une quelconque structure sauf que les clusters tunisiens de plus de deux individus se retrouvent bien ensemble (fig. 32) avec un cluster du Tessin (Cen16 qui comporte d’ailleurs deux mâles et une femelle). Ceci ne contredit pas que l’effet Wahlund pourrait être issu de la présence dans chaque site d’individus issus des mêmes pontes. Ceci implique une forte variance du succès de survie entre pontes : de nombreux individus issus seulement de quelques pontes accèdent à l’âge adulte (voir Chevillon et al., 2007a, pour un résultat similaire sur la tique du bétail).

Figure 32

Commentaires sur l’analyse des fichiers par BAPS

78Contrairement à ce qui pourrait être suggéré à la lecture du manuel d’utilisation de BAPS, les clusters obtenus ne présentent pas ici une structure panmictique, mais conservent un déficit important d’hétérozygotes sauf pour IR08. Nous verrons, avec les analyses suivantes, que ces clusters reflètent probablement en grande partie une réalité biologique de nature assez complexe (races d’hôte, structures familiales), et qui devra conduire à d’autres études. BAPS ne représente ici qu’un outil d’argumentation et d’orientation de futures investigations, pas un générateur de vérités.

Conclusion sur les déficits en hétérozygotes

79À l’occasion de ces premières analyses, nous pouvons constater qu’une analyse d’un jeu de données de génétique des populations requiert de la patience, de la méthode, ainsi qu’une bonne batterie de tests. Il était cependant nécessaire d’aller jusqu’au bout avant d’aller plus loin. Nous savons maintenant que ces tiques sont structurées à une échelle locale, ce qui explique une grande partie des déficits en hétérozygotes. Cet effet Wahlund résulte probablement d’une structure en groupes familiaux. L’existence d’espèces cryptiques n’est en effet pas soutenue par nos analyses NJTREE ni par l’absence totale de déséquilibre de liaison. Nous savons également qu’une partie non négligeable de ces déficits provient de l’existence d’allèles nuls (Loci IR25, IR32, IR39). Pour ces derniers, il y a donc un risque de surestimer la différenciation entre sous-échantillons, mais seulement pour des niveaux de différenciation atteignant au moins 10 % (F_ST = 0,1), en dessous de quoi l’effet devient faible (Chapuis et Estoup, 2007). Nous verrons que les niveaux de différenciation entre populations d’I. ricinus se trouvent bien en dessous de cette frontière. Enfin, un locus (IR27) a montré des évidences de dominance des allèles courts. Dans la mesure où ce phénomène modifie l’hétérozygotie et l’estimation des fréquences alléliques, il faudra être constamment vigilant quant aux résultats obtenus par la suite. Si nous avions un locus de plus sans allèle nul, j’aurais même conseillé de le supprimer. Ce n’est malheureusement pas le cas. Il faudra juste vérifier que chaque résultat ultérieur n’est pas sous la dépendance de ce seul locus. L’idéal aurait été d’avoir sept loci comme IR08, mais non liés à l’X ! Mais on ne choisit pas et les problèmes de marqueurs chez les parasites et vecteurs représentent un souci récurrent.

80Une autre conclusion importante est qu’un déficit en hétérozygotes non entièrement expliqué par des allèles nuls exclut les causes endogamiques (croisements frère/sœur, autofécondation…) qui tendent à augmenter l’homozygotie et donc à dévoiler les homozygotes nuls (blancs), d’une part, et suggère, d’autre part, plutôt un effet Wahlund, qui augmente la diversité génétique H_s sans augmenter l’hétérozygotie observée (d’où augmentation du F_IS, cf. équation 19 en première partie de ce manuel, p. 49). Dans le cas d’un effet Wahlund, il est donc normal que les procédures de détection d’allèles nuls ne suffisent pas à expliquer entièrement les déficits en hétérozygotes, même si ces derniers sont présents, comme l’attestent la présence fréquente d’individus blancs, ainsi que la forte variance du F_IS entre loci.

Structure des populations et schémas de différenciation

81Nous avons ici une espèce à sexes séparés. La première chose à tester est s’il n’existe pas une différence entre femelles et mâles tiques, liée par exemple à un biais de dispersion spécifique de chaque sexe (Goudet et al., 2002 ; Prugnolle et De Meeûs, 2002). En plus, on sait que c’est probablement le cas ici puisque ce signal fut détecté précédemment (De Meeûs et al., 2002a), mais aussi lors de notre recherche de dominance des allèles courts. Par ailleurs, il est intéressant de vérifier si le signal persiste en tenant compte de la microstructure en clusters, même s’il a été montré que celle-ci a peu (pas) d’effet sur la structure à plus large échelle, si la microstructure n’est pas trop forte (Fontanillas et al., 2004).

Structure génétique spécifique à chaque sexe des données brutes (sans tenir compte de BAPS)

82Comment suspecter qu’un biais de structuration existe entre mâles et femelles ? Soit en effectuant directement le test “Sex biased dispersal” de Fstat, soit, comme cela a été le cas pour les données présentes, en testant la différenciation locale entre tiques mâles et femelles. La justification de ce test est qu’un tel signal avait été suggéré chez cette espèce en Irlande pour un locus enzymatique (Healy, 1979). Nous allons donc mesurer et tester la différenciation entre mâles et femelles de chaque sous-échantillon. Pour ce faire, il faut construire un fichier Fstat (par exemple) où femelles et mâles de chaque site sont considérés comme appartenant à des échantillons différents. Appelons ce fichier “IRTotBrutSexBias.dat” et chargeons-le dans Fstat (après avoir ouvert Fstat il faut aller dans “File”, “Open”, etc.). On coche ensuite dans Fstat l’option “Fst per pair of samples” et la case “Pairwise tests of differentiation”, ainsi que la case “5/100” du “Nominal level for multiple tests”. Cette dernière case est choisie pour gagner du temps. Ici, Fstat donnera le seuil Bonferroni de significativité corrigé par le nombre de tests (276 ici). À ce seuil, une P-value sera significative si elle est inférieure ou égale à 0,05/276 = 0,00018 et Fstat ajuste le nombre de permutations nécessaires pour atteindre cette valeur, soit 5 520, ce qui est bien suffisant. Avec “1/100” on obtient 27 600, ce qui est beaucoup. En plus, à ce niveau, le Bonferroni est beaucoup trop conservateur. De toutes façons, comme nous n’allons utiliser que les résultats par paire locale de femelles et de mâles, nous n’appliquerons pas cette procédure. Après avoir cliqué sur “Run” et attendu la fin des permutations, deux fichiers sont à consulter. “IRTotBrutSexBias.fst” donne les F_ST par paire et “IRTotBrutSexBias-pp.pvl” donne les P-value du test de randomisation des génotypes par paire de sous-échantillons. Dans ces fichiers, il faut garder les valeurs correspondant aux paires femelle-mâle de chaque site-année. Si vous faites cela, deux probabilités sont significatives sur les 12 tests (17 %). Un test binomial peut alors être effectué sous R avec la commande suivante :

83binom.test(2, 12, p=0.05, alternative="greater")

84Le test est unilatéral, car on regarde si le nombre 2 n’est pas plus grand qu’attendu sous l’hypothèse nulle. Le test donne une P-value de 0,12, ce qui n’est pas vraiment significatif, mais témoigne d’un signal possible. Le test généralisé entrepris avec MultiTest et un k’ = 6 donne une P-value globale seuil de 0,6015 (La notice d’utilisation de ce programme est suffisamment détaillée pour ne pas avoir à reproduire ici un tuteurage pas à pas).

85Pour effectuer le véritable test de structuration sexe-spécifique, il faut remanier quelque peu le fichier initial des données afin de le mettre au format requis par Fstat pour l’analyse du biais de dispersion sexe-spécifique (Biased dispersal menu). Référez-vous à l’aide en ligne de Fstat pour construire ce fichier. Nous allons nous focaliser sur les échantillons 1996 uniquement. Une fois ce fichier constitué, il faut lancer Fstat, cliquer sur le menu “Biased dispersal” et y charger le fichier requis. Il faut ensuite sélectionner le test “Two sided” (on n’a en principe pas de préjugé pour l’instant) et cocher toutes les options comme dans la figure 33.

Figure 33

86Vous remarquerez dans la figure 33 que les cases du F_IS et du H_o sont cochées comme les autres, alors que cela n’a aucun sens. En effet, puisque nous avons codé les mâles homozygotes au locus IR08, il y aura nécessairement une différence mâle femelle à ce niveau. Cependant, quand cette option n’est pas cochée, on perd une partie de l’information sur H_s dans le fichier de sortie. Il conviendra donc, dans ce fichier, d’ignorer les résultats sur H_o et F_IS. Le logiciel crée cinq fichiers, trois fichiers .dat au format Fstat (les données totales, les femelles, les mâles), le fichier de permutations et le fichier .res des résultats (le plus utile). Ces derniers indiquent que les femelles sont bien mieux assignées que les mâles (IA_c = 0,36 et IA_c = – 0,56 pour les femelles et les mâles respectivement, P-value = 0,0005) et que les femelles sont localement mois diverses génétiquement (H_s = 0,79) que les mâles (H_s = 0,81) (P-value = 0,027), ce qui va dans le sens d’un biais de dispersion mâle (les femelles disperseraient moins). Par contre, le F_ST et la variance d’assignement répondent en sens inverse (mais non significativement heureusement). Pourtant, ce sont ces derniers paramètres (F_ST et variance d’assignement) qui doivent théoriquement signaler les premiers un biais de dispersion (qui donnent les tests les plus puissants) (Goudet et al., 2002). Nous discuterons de ce paradoxe plus loin.

87Afin de tester si la Tunisie n’est pas responsable seule de ce résultat, recommençons avec les données de Suisse 1996. Dans ce cas, on a des résultats comparables avec une P-value = 0,0004 pour l’assignement, mais une P-value = 0,06 marginalement significative pour H_s. Cantonnons-nous (normal pour la Suisse) au Plateau Suisse en excluant le site Monte-Ceneri du Tessin. Cette fois, les P-values tombent à 0,0002 et 0,02 pour les assignements et H_s respectivement. En restreignant l’échantillonnage aux sites du nord-ouest de la Suisse (il faut supprimer les sites Gorges-du-Trient et Dorénaz), sans oublier de le signaler en en-tête du fichier de données (il n’y a plus que cinq sites), on obtient une confirmation de ce qui était observé (tabl. 11), mais sur une échelle plus réaliste quant aux interprétations biologiques (en fin de ce chapitre). Il semble donc bien y avoir un biais de dispersion mâle (ou à tout le moins un biais de structuration génétique en faveur de ces femelles). En retirant chaque locus un à un et en recommençant l’analyse (donc cinq traitements), vous pourrez vérifier qu’aucun locus n’est responsable à lui seul du signal. On constate même, pour les données sans IR08, que le F_IS est significativement supérieur chez les mâles (tabl. 11). On pourra ici se contenter de refaire ces analyse sur les échantillons du Nord-Ouest et en unilatéral pour compenser la perte de puissance. La question qui se pose ensuite est de savoir si tenir compte des résultats de BAPS (microstructuration) change cette conclusion. Pour ce faire, il faut réanalyser les données en tenant compte des clusters définis par BAPS.

Structure génétique spécifique à chaque sexe des données clusterisées par BAPS

88Nous prendrons ici le fichier de données 1996 de Suisse uniquement et les clusters obtenus en ne séparant pas les mâles des femelles (évidemment). Il faudra prendre garde à ne garder que les clusters contenant au moins une femelle et un mâle, car sinon Fstat va planter (comme on dit). Nous allons dans un premier temps effectuer l’analyse sur tous les clusters de tous les sites. Le label “Pop” va donc se positionner entre chaque cluster. On peut faire le test en unilatéral, mais au vu des résultats vous verrez vite qu’il convient de repartir sur une base de tests bilatéraux. Les résultats sont en effet spectaculairement divergents des précédents (tabl. 12).

89Ce résultat, très déconcertant au premier abord, est sous très forte influence du locus IR08, bien que les autres loci répondent dans le même sens (sauf peut-être IR32). Comme il s’agit peut-être d’un phénomène local, nous allons refaire les mêmes analyses, mais dans chaque site de 1996 (y compris la Tunisie). Le résultat des tests sur le F_ST figure dans le tableau 13. Le signal reste le même, mais semble disparaître sans le locus IR08. Il se pourrait que ce locus soit diagnostique de certains groupes de tiques. Pour vérifier cela, il faut reprendre le fichier initial de données et grouper les individus, dans chaque site, selon leur génotype au locus IR08. Ce faisant, on recalcule sur cette nouvelle partition le F_IS et le F_ST avec Genetix, ce qui donne 0,47 et 0,02 respectivement, alors qu’on attend un faible F_IS et un fort F_ST. IR08 n’est manifestement diagnostique de rien du tout et le fait qu’il donne les meilleurs résultats provient vraisemblablement de sa qualité (peu ou pas d’allèles nuls et très faible variance des différents estimateurs).

90Il y a donc manifestement un effet cluster que nous essayerons d’interpréter plus loin. Afin de vérifier quand même si notre biais de dispersion spécifique femelle existe toujours même en tenant compte de l’effet Wahlund présent au sein de chaque site, la solution qui nous reste consiste à ne garder qu’un seul représentant ou une femelle et un mâle par cluster dans chaque site (nord-ouest de la Suisse 1996). On prendra le premier des individus ayant le génotype le plus complet de chaque cluster afin de conserver le plus de puissance possible. Par exemple, si dans un cluster d’un site quelconque, il n’y a que des mâles on ne prend qu’un individu, si possible génotypé aux cinq loci. Même chose pour des clusters de femelles. Pour les clusters mixtes, on prend la première femelle la plus complète et le premier mâle le plus complet. On obtient ainsi un jeu de données de cinq sites avec un nombre d’individus fortement réduit par site. C’est aussi la raison pour laquelle les tests seront unilatéraux (les femelles dispersent moins). Le résultat de cette analyse figure dans le tableau 14 où on retrouve bien le signal initial suggérant un biais de dispersion femelle, à la différence que tous les paramètres vont dans le bon sens, même si c’est toujours AI_c qui donne la seule P-value significative.

Interpréter l’ensemble des résultats sur les biais de structuration

91Il semble bien y avoir un biais de dispersion biaisé pour les femelles (elles disperseraient moins) à l’échelle du plateau Suisse (ou même de régions plus restreintes), mais le signal est brouillé par une micro-structuration qui existe localement. Le fait que dans chaque site, les clusters trouvés par BAPS contiennent des femelles beaucoup plus hétérogènes que les mâles à l’intérieur de chaque cluster, mais beaucoup moins différentes d’un cluster à l’autre peut être interprété de deux façons. La première suggérerait que le biais de dispersion spécifique à chaque sexe s’inverse à petite échelle, mais on ne voit pas bien comment. La seconde suppose que les clusters correspondent plus ou moins à des frères et sœurs issus d’une même ponte et que les femelles ont une réussite beaucoup plus homogène que les mâles. Ne parviendraient à l’âge adulte, selon cette hypothèse, que beaucoup de mâles par ponte, mais de peu de pontes, alors que les femelles représenteraient un échantillon plus aléatoire des pontes (moins de sœurs que de frères dans chaque site). Pour confirmer cette interprétation, une approche théorique de modélisation/simulation serait nécessaire, mais dépasserait alors le cadre ambitionné par cet ouvrage. Enfin, ces clusters pourraient correspondre à des cohortes différentes (chevauchement de générations), très différenciées (dérive forte) et cela surtout pour les mâles dont beaucoup viennent d’ailleurs. Ici aussi, une approche théorique s’avérerait nécessaire. Il est cependant raisonnable d’imaginer que si les larves et les nymphes mâles sont plus souvent retrouvées sur des hôtes très dispersants, alors il y a de fortes chances que chacun de ces individus hôtes porte des mâles apparentés (surtout les larves). Une fois dispersé et gorgé, chaque groupe a une chance très inégale de trouver un habitat favorable à la mue suivante. Il en résulterait que seuls certains groupes, parfois composés d’individus très apparentés (frères), survivraient dans une zone éloignée de leur site d’éclosion, alors que beaucoup de groupes mâles seraient éliminés. Si les larves et nymphes femelles préfèrent, quant à elles, les hôtes peu dispersants (petits rongeurs), il est probable que la survie de ces femelles soit distribuée plus aléatoirement entre femelles de pontes différentes. Ceci pourrait au final expliquer notre effet Wahlund produit en majorité par les tiques mâles.

Différenciation globale et isolement par la distance

92Plusieurs éléments nous incitent ici à manquer d’optimisme. Il y a en effet de nombreux allèles nuls, un effet Wahlund local, de la dominance d’allèles courts à un locus, sans parler d’autres problèmes mis en évidence lors d’études de pedigrees (De Meeûs et al., 2004a). Si on ajoute à cela que manifestement un biais de dispersion spécifique à chaque sexe existe, supposant qu’un des deux sexes migre beaucoup (voir Goudet et al., 2002) et donc qu’une faible structuration en résulte nécessairement, la probabilité de trouver une structuration génétique devient faible, et c’est un euphémisme. Nous allons quand même tenter notre chance, et ce pour plusieurs raisons. D’abord, parce que nous ne sommes pas arrivés jusqu’ici pour se mettre à bailler aux corneilles, ensuite parce que « c’est la nuit qu’il est beau de croire à la lumière » (Rostand, 1908).

Définir différents niveaux de subdivision pour l’analyse hiérarchique

93Nous ne considérerons ici que les échantillons de 1996. Nous pouvons envisager, grâce à HierFstat (Goudet, 2005), n’importe quelle structure du moment que cette dernière reste hiérarchique. Nous allons donc dans un premier temps considérer (référez-vous au besoin à la figure 15) l’Europe-Afrique comme tout, suivi de la Tunisie versus la Suisse, puis le Tessin versus le nord des Alpes et enfin le groupe Gorges-du-Trient, Dorénaz contre le plateau Suisse (Eclepens, Montmollin, Neuchâtel, Staadswald, Bern). Référez-vous à De Meeûs et Goudet (2007) pour des détails sur la confection d’un fichier HierFstat.

Analyse hiérarchique sur données brutes (pas de cluster BAPS)

94Il faut donc créer un fichier avec quatre (hiérarchie) plus cinq (loci) colonnes. La première colonne correspond donc au continent, Cont avec 1 l’Europe (= la Suisse, et alors ?) et 2 pour l’Afrique (Tunisie). La deuxième colonne (NrdWTessin) va coder pour l’appartenance aux cantons du nord et nord-ouest de la Suisse (1), pour celle du Tessin (2) (Monte-Ceneri) ou la Tunisie (3) qui n’est pas plus subdivisée, mais doit être aussi codée dans cette colonne. La troisième colonne (NrdWNS) correspond à l’appartenance ou non au nord-ouest (1) ou au sud-ouest (Gorges-du-Trient, Dorénaz = 2) de la zone du nord des Alpes suisses. Le Tessin et la Tunisie étant codés 3 et 4 respectivement dans cette colonne. La quatrième colonne (Site) correspond aux sites eux-mêmes (1 à 9). Les cinquième à neuvième colonnes correspondent aux cinq loci, le premier, IR08, étant codé homozygote pour les tiques mâles. Appelons le fichier ainsi construit “IRTot96HierFstat.txt”. L’analyse va se faire sous HierFstat 0.04-4 (Goudet, 2006, mis à jour de Goudet, 2005) comme décrit dans De Meeûs et Goudet (2007). N’oubliez pas de remplacer les données manquantes “000000” par “NA”. Lancez le logiciel R. Chargez le package HierFstat (Menu “Package”, “Chargez le package”, “hierfstat”). Changez de répertoire pour travailler dans celui où le fichier de données “IRTot96HierFstat.txt” se trouve (Menu “Fichier”, “Changer le répertoire courant”). Dans la console R, tapez la succession de commandes (chaque ligne correspond à une commande devant être suivie d’un retour charriot), en respectant les majuscules et minuscules (distinctes en langage R):

95Cette dernière commande produit le résultat suivant :

96Dont l’interprétation est la suivante :

97F_IS = 0,459 (nous retrouvons ici un résultat ancien et sans valeur, car les mâles sont artificiellement homozygotes ici au locus IR08), F_Site/NrdWNS = 0,0028, F_{NrdWNS/NrdWTessin} = 0,0025, F_{NrdWTessin/Cont} = 0,0057 et F_Cont/Total = 0,106. Toutes ces valeurs de différenciation sont très faibles sauf pour la Suisse contre la Tunisie. Il faut tester ensuite la significativité de ces différentes partitions en commençant par la plus incluse, le site :

98On voit bien que le site (comme on le craignait) n’influence en rien la partition de l’information génétique. Nous allons donc supprimer ce facteur de la hiérarchie :

99Le facteur NrdWNS, séparant les sites Dorénaz-Gorges-du-Trient de l’ensemble des sites suisses du Nord-Ouest, ne semble pas influencer davantage la structure génétique des tiques. Si nous le supprimons à son tour, nous obtenons :

100Si nous choisissons de garder le facteur NrdWTessin (marginalement significatif, P-value = 0,058) cela aboutit à :

101Si on élimine le facteur NrdWTessin, il faut alors repasser par Fstat. Il n’y a en effet plus que trois niveaux hiérarchiques avec deux sous-populations représentées par l’ensemble des tiques suisses, d’une part et par celles de Tunisie, d’autre part. On aboutit à un F_ST = 0,113 très significatif (P-value < 0,0001) entre les tiques de Suisse réunies en une seule population et la Tunisie.

102Avec un H_s = 0,832, cela correspond à un F_ST’ = F_ST/F_STmax = 0,113/(1 – 0,832) = 0,673, ce qui est relativement considérable et témoigne du peu de migration entre les deux pays. Par contre, à l’échelle de la Suisse, cette migration est forte et même si les Alpes apparaissent comme un facteur limitant, tout semble se passer comme si, génétiquement au moins, on avait à faire à une seule unité à cette échelle.

103Qu’en est-il si nous tenons compte des clusters trouvés par BAPS ?

Analyse hiérarchique sur données clusterisées par BAPS

104Nous allons donc utiliser le fichier de données précédent avec une colonne supplémentaire correspondant aux clusters trouvés avec BAPS. En suivant alors une procédure identique à celle décrite plus haut, nous pouvons constater que les facteurs ClusterBAPS (F_Clust/Site = 0,3, P-value = 0,001, ce qui, il faut bien l’avouer, est trivial) qui mesurent la partition génétique entre clusters d’un même site, et Continent (F_{Continent/Total} = 0,11, P-value = 0,001) qui mesure la différenciation entre Suisse et Tunisie, constituent les deux seuls facteurs qui structurent les sous-échantillons de façon significative.

105Si nous ne gardons qu’un mâle ou une femelle ou un individu par cluster, comme pour le tableau 14, le résultat de l’analyse par HierFstat ne montre plus aucune différenciation, à moins d’ignorer tous les facteurs sauf le continent (analyse par Fstat, F_ST = 0,09, P-value = 0,001).

Test d’isolement par la distance

106Nous ne travaillerons ici que sur les échantillons suisses de 1996. D’abord parce que la Tunisie est trop éloignée par rapport aux distances entre échantillons suisses. Il y aurait deux groupes de points. Procéder à un test de régression entre deux points n’a pas de sens, le plus court chemin entre eux étant nécessairement une droite, c’est dans tous les bons livres de statistiques. Or, le test d’isolement par la distance est une forme de régression où on cherche à expliquer une différence génétique croissante par un éloignement géographique. Ensuite, il n’y a pas assez d’échantillons en 1995.

107Pour le test, il faut configurer un fichier avec deux demi-matrices, l’une pour les distances géographiques entre paire de sites et l’autre pour les F_ST (estimés par θ) correspondants. Pour les distances géographiques, vous pouvez vous aider de la figure 15. Pour les F_ST, il suffit de prendre la sortie “IRTot96CH.fst” que Fstat a produit en analysant le fichier “IR96CH.dat” des données suisses 1996, si vous avez toutefois coché la case “Fst per pair of samples”. En ce qui me concerne, j’obtiens les matrices représentées dans le tableau 15. Le test va être effectué selon la méthode décrite par Rousset (1997) pour un schéma en deux dimensions. Nous allons donc effectuer un test de Mantel sur la corrélation entre le F_ST/(1 – F_ST) et le log népérien (ou naturel) de la distance géographique. Nous allons utiliser Genepop 3 pour faire ce test et donc formater les données dans ce sens et les sauvegarder dans un fichier que nous appellerons IR96CH.mig. Ce fichier doit être configuré comme présenté dans la figure 34.

Figure 34

108L’étape suivante consiste à lancer Genepop 3. Éviter de double cliquer sur le fichier genepop.bat, mais préférez ouvrir une session DOS en lançant une “Invite de commandes” dans le menu “Accessoires” de Windows. Dans la fenêtre DOS, et si Genepop est dans le répertoire “Genepop” du disque D, tapez “D:”, puis “Entrée”, puis “cd Genepop”, puis “Entrée”. Vous êtes dans le répertoire Genepop. Tapez alors “isolde”, puis “Entrée” pour lancer le programme d’isolement par la distance. À l’invite, tapez le nom complet du fichier de données puis “Entrée”. Le logiciel vous demande ensuite quel type de distance (non transformée ou Log) et quel type de mesure de différenciation vous souhaitez tester (X, qui figure dans la matrice ou X/(1 – X) ). À vous de choisir la méthode appropriée. Genepop vous demande ensuite la distance minimale en deçà de laquelle la mesure de corrélation ne tient plus compte des données, car en deçà d’un certain niveau la réponse a en effet tendance à ne plus suivre un modèle clair (Rousset, 1997). Réfléchissez à ce que devrait être cette distance minimale et tapez-la. Le nombre de randomisations vous est ensuite demandé. Tapez 1 000 000 pour être sûr d’obtenir une P-value suffisamment précise. Genepop vous demande, comme pour tous ses tests passant par randomisation, quatre nombres qui serviront de « graines » pour générer des nombres pseudo-aléatoires qui conditionnent le départ des randomisations. Tapez ce que vous voulez entre 1 et 168, comme indiqué avec un retour charriot après chaque chiffre. Quand les randomisations sont terminées, Genepop le signale avec un bip assez désagréable, mais qui ne doit pas vous effrayer (en général, je coupe le son avant). D’après une collègue avisée (TG), il n’y a pas de son sur la version Web du programme qui, par ailleurs, est sensiblement plus conviviale. Genepop a alors créé deux fichiers, l’un porte l’extension .ISO qui vous donne le résultat du test avec les paramètres de la régression et la P-value. Le second porte l’extension .GRA et donne les coordonnées en colonne de cette régression afin de pouvoir dessiner le graphique correspondant, comme représenté dans la figure 35. On y remarque que la relation n’est pas significative. Il semble cependant qu’une tendance existe. Peut-être l’existence d’une sous-structure nuit-elle à la clarté du signal ? Pour le vérifier, il suffit de procéder au même test, mais avec les données où un ou deux (de chaque sexe) individus par cluster avaient été gardés (voir p. 156-158). Le résultat change du tout au tout puisque la relation devient très significative, comme en témoigne la figure 36. Ceci permet de calculer le voisinage Nb = 1/b = 173 individus (Watts et al., 2007), le produit de la densité d’individus par km² par la surface de dispersion des descendants reproducteurs par rapport à leurs géniteurs, en utilisant la méthode de Rousset (1997) décrite en première partie (p. 90-92), ce qui donne Dσ² = 1/(4π0,00577) = 13,78. On peut aussi estimer le nombre d’immigrants présents dans une sous-population (Rousset, 1997), Nm = 1/2πb = 28 individus par génération. Il serait dommage de s’arrêter ici et nous allons donc essayer d’appréhender, même grossièrement, dans quelle gamme de valeurs se situe la densité de tiques afin d’en pouvoir extraire la surface de dispersion de ces tiques.

Figure 35

Figure 36

Estimation d’effectifs efficaces, extrapolation des densités et de la dispersion

Effectifs efficaces des tiques de Suisse

109Pour des raisons de commodité et de cohérence (les tiques tunisiennes n’ayant pas été échantillonnées de la même manière), nous nous focaliserons sur les échantillons de Suisse. Rappelons-nous que nous avons rencontré de gros déficits en hétérozygotes (allèles nuls et dominance d’allèles courts), ainsi que la présence d’un fort effet Wahlund. Nous ne travaillerons donc qu’à partir de méthodes indépendantes de l’hétérozygotie, telle que celle proposée par Bartley et al.(1992), basée sur les déséquilibres de liaison et implémentée par NeEstimator (Peel et al., 2004). Nous n’utiliserons que les données clusterisées par BAPS où seuls subsistent un ou deux individus par cluster dans chacun des huit sites suisses pour éviter l’effet confondant dû à l’effet Wahlund.

110Il faut créer un fichier par site dans un format proche de Genepop comme dans la figure 37.

Figure 37

111Il faut ensuite lancer le programme NeEstimator (après l’avoir installé sur votre machine, bien entendu). Une fenêtre d’avertissement sur le copyright et sur la manière idoine de citer ce logiciel apparaît. Cliquez sur OK pour accéder au programme qui apparaît dans une fenêtre comme dans la figure 38. Comme indiqué sur la figure 38, cliquez sur le menu déroulant “File” et “Open”, ce qui permet d’ouvrir la fenêtre “Analysis”.

Figure 38

112Dans la fenêtre “Analysis”, une série d’onglets apparaît et vous positionne sur celui du format de vos données “Data Format” où il n’y a rien à changer, car vous avez choisi le format par défaut. Allez à l’onglet “Data Files”. Là il n’y a qu’un seul bouton “Load” qui vous permet de charger votre jeux de données, ce que vous faites (fig. 39). Une fois que vous avez choisi le fichier, le logiciel vous demande à quelle génération ces données correspondent-elles. Laissez la valeur par défaut “0”, car nous n’utiliserons pas ici la méthode des moments de Waples (1989) (cf. p. 107 en première partie) et cliquez sur “OK”. Dans le menu déroulant “NeEstimator”, cliquez sur “Run” (fig. 40). Ce qui fait apparaître un message qui vous avertit qu’avec un seul échantillon, on ne peut utiliser les méthodes temporelles “Moment based” et vous demande si vous souhaitez continuer avec les méthodes à un seul échantillon. Vous répondez “Oui” bien entendu. Le résultat est affiché sous forme de tableau que je vous conseille de sauvegarder au format NeEstimator (NeA). Je conseille aussi de transcrire tous les résultats dans un tableur au fur et à mesure afin de disposer de l’ensemble dans un seul fichier. C’est ce qui est représenté dans le tableau 16.

Figure 39

Figure 40

113Ici, bien que nous disposions d’échantillons espacés dans le temps (Bern, Gorges-du-Trient et Staadswald), ces échantillons ne sont séparés que d’une année, soit environ 1/3 du temps de génération d’I. ricinus. Ici, les adultes présents d’une année sur l’autre font partie de cohortes séparées et qui, même à long terme, auront du mal à échanger des gènes. La différenciation entre ces cohortes, déjà remarquée par De Meeûs et al. (2002a), va tendre à être très supérieure à celle qui existe réellement entre deux générations d’adultes reproducteurs. L’utilisation des méthodes temporelles sur nos données aboutira donc à de fortes sous-estimations des effectifs efficaces. Faites-le et vérifiez qu’effectivement, compte tenu qu’il n’y a qu’un tiers de génération séparant 1996 de 1995, les estimations obtenues par la méthode de Waples donnent des effectifs efficaces proches de 0, ce qui n’est pas très conforme à la perception que l’on peut avoir sur le terrain.

114En reprenant le tableau 16, nous obtenons par conséquent un effectif efficace de 268 en moyenne sur l’ensemble des échantillons avec un intervalle de confiance à 95 % de [119, 541], avec des valeurs minimales et maximales de 73 et 700 respectivement. Ces nombres paraissent plausibles, compte tenu de l’effet Wahlund reflétant probablement un fonctionnement particulier des populations de tiques susceptible d’en réduire sensiblement l’estimation de leurs effectifs efficaces.

115En reprenant les données avec un ou deux individus par cluster BAPS, les valeurs obtenues sont plus grandes en moyenne (596) avec un minimum et un maximum de 75 et 1 057 respectivement3.

Extrapolation des densités et des distances de dispersion des tiques en Suisse

116Il faut dans un premier temps estimer sur quelle surface se distribuent les tiques. Ici, c’est difficile et on ne peut pas dire grand-chose de plus que les surfaces d’échantillonnage s’étendaient grossièrement sur S = 0,2 km². Ceci signifie (mais vous vous en doutiez probablement) que les estimations à venir seront tout à fait approximatives. À partir de là, les densités sont faciles à calculer (N_e/S). La densité moyenne devient 1 340 tiques reproductrices/km² 95 % CI = [594, 2 706] avec un minimum et un maximum de 367 et 3 502 tiques/km² respectivement (tabl. 16). En réutilisant les résultats de la régression de l’isolement par la distance D_eσ² = 13,78 (voir p. 166), on aboutit à une surface de dispersion moyenne entre adultes et leurs parents d’environ 0,01 km² [0,005, 0,023] avec un minimum et un maximum de 0,004 et 0,038 km² respectivement. Autrement dit, la distance moyenne séparant un adulte reproducteur de ses géniteurs est d’un ordre de grandeur de 100 m par génération (donc tous les trois ans environ), un intervalle de confiance à 95 % de bootstrap = [71, 152] et un maximum et un minimum de 63 à 195 m, ce qui est relativement modeste. Les données clusterisées par BAPS conduisent à une densité de 3 000 tiques par km² et une dispersion de moins de 60 m par génération. Donc, sachant que l’estimateur sans doute le moins biaisé est le produit D_eσ², la dispersion par génération est, quoi qu’il en soit, extrêmement modeste à moins d’évoquer des densités (effectifs) efficaces extrêmement faibles. Il en va donc de même en ce qui concerne la propagation des maladies par les tiques.

Conclusions de la 1^re édition de ce manuel sur la biologie et la génétique des populations d’I. ricinus en Suisse

117Il existe un déficit important en hétérozygotes dans les populations d’I. ricinus (F_IS = 0,39) dont une majeure partie (64 %) est expliquée par un effet Wahlund important.

118Le F_IS = 0,14 résiduel correspondrait à du « stuttering », à de la dominance d’allèles courts et à des allèles nuls. Pour tester les allèles nuls dans les clusters de BAPS, on ne peut pas utiliser Micro-Checker (échantillons trop petits). Nous pouvons néanmoins tester s’il existe une relation positive entre le nombre de blancs à un locus et le F_IS à ce locus. En effet, en reprenant les données clusterisées et en séparant les mâles des femelles en deux fichiers, il est facile de compter les blancs pour chaque locus avec la fonction “SI” d’Excel. Il suffit de créer autant de nouvelles colonnes qu’il y a de loci et de remplir chacune avec les instructions de type “= SI(G2 = “000000”) ; 1;0)” pour inscrire “1” quand on a un blanc. À la fin de chacune de ces colonnes, on tape une instruction du type “= somme(L2:L147)” pour obtenir la totalité des blancs à ce locus sur l’ensemble des clusters. Le F_IS de chaque locus est récupérable dans les deux fichiers de sortie Fstat de l’analyse des deux jeux de données (un pour le femelles et un pour les mâles) avec les données clusterisées par BAPS, que j’ai personnellement nommés IRTotBAPSClustMalManqIR08Females.dat et IRTotBAPSClustMalManq IR08FMales.dat respectivement, et où on aura pris soin d’éliminer le locus IR08 du fichier des mâles. Quand on a fait ceci pour les femelles et les mâles, on obtient le jeu de données présenté dans le tableau 17. La corrélation entre le nombre de blancs et le F_IS peut être analysée par un test de corrélation de Spearman (test non paramétrique). Ce test est facile à réaliser sous R. Si le fichier de données correspondant au tableau 17 s’appelle “AllelesNulsClustersBAPS.txt”, alors il suffit de lancer R, et de se placer dans le répertoire contenant ce fichier (menu déroulant “Fichier”, “Changer le répertoire courant”).

119Ensuite, il faut taper les instructions suivantes :

120ce qui renvoie au résultat :

121La figure 41 illustre la relation positive forte entre les deux variables.

Figure 41

122Nous pouvons également tester de nouveau la dominance d’allèles courts au locus IR27 en prenant les F_IS par allèle donnés par la sortie Fstat sur les mâles et les femelles séparément. Nous ne distinguerons en effet pas les clusters, car il y en a trop. Sous R, la procédure est comparable à celle utilisée en p. 135-140, sauf que nous n’utiliserons que le polynôme d’ordre deux de la taille des allèles et le sexe comme variables explicatives du F_IS. Le résultat est de nouveau très significatif et on explique donc toujours une forte proportion du F_IS par ce phénomène de dominance des allèles courts, comme illustré par la figure 42. Suivent les instructions R :

123ce qui renvoie au résultat suivant :

Figure 42

124Comme cela a été vu au début de ce paragraphe, la majeure partie (64 %) du F_IS initial est expliquée par un effet Wahlund. Cet effet Wahlund est plus prononcé chez les mâles dont les clusters se trouvent plus différents entre eux que les femelles. Dans l’hypothèse de spécificités d’hôtes différentes des larves et/ou nymphes femelles et mâles, des groupes de larves ou nymphes mâles fortement apparentés seraient transportés ensemble sur le même hôte très dispersant (oiseau) avec de forts risques de tomber sur des sites défavorables lors du détachement, à la fin du repas sanguin. Les mâles retrouvés adultes dans nos échantillons correspondraient alors aux quelques groupes d’apparentés ayant eu la chance de tomber ensemble dans un site favorable. Les larves ou nymphes femelles seraient, quant à elles, plus souvent retrouvées sur des hôtes très peu dispersants, comme des petits rongeurs très territoriaux. Il en résulterait un apparentement réparti beaucoup plus aléatoirement pour les femelles dans chaque site. Il y a un fort biais de dispersion spécifique à chaque sexe (les femelles dispersent très peu). Ce biais est partiellement masqué par l’effet Wahlund, et il est plus facilement visible quand cet effet est corrigé (données réduites), et l’indice d’assignement corrigé AI_c semble à cet égard beaucoup plus robuste que sa variance vAI_c et le F_ST.

125Cet effet Wahlund nuit considérablement à l’image perçue au niveau de la structuration à l’échelle de la Suisse. Quand cet effet est contrôlé (au moins en grande partie), on observe un isolement par la distance très significatif, et les adultes non gorgés d’I. ricinus paraissent distribués en populations locales de tailles importantes (plus de 1 000 tiques par km²) et se dispersant difficilement à plus de 200 m par génération.

126Il reste cependant bien d’autres questions et toutes ces hypothèses doivent être testées sur le terrain. Cette étude ouvre de nombreuses et prometteuses perspectives de recherche que je vous laisse le soin de discuter.

Discussion des résultats obtenus par des méthodes plus récentes d’analyse pour la 2^e édition de ce manuel

127Pour cette réédition j’ai préféré adopter une approche quelque peu différente. J’ai choisi en effet de combiner les résultats obtenus avec deux fichiers, l’un sans IR08 mais avec tous les individus et l’autre sans les mâles mais avec tous les loci. Je n’ai pas partitionné les sous-échantillons en fonction du sexe des tiques. J’ai compilé l’ensemble des résultats dans le fichier « IxodesResults.xlsx » que vous pouvez télécharger sur mon site web :  http://www.t-de-meeus.fr/Data/DataLivreInitiation/Data.html.

128On y remarque que de nombreuses conclusions restent inchangées, mais que d’autres sont à revoir. Il semble que l’hypothèse d’un effet Wahlund n’apparaît pas expliquer grandement les données. Je le soupçonnais déjà, car un effet Wahlund important aurait dû générer davantage de déséquilibres de liaisons dans cet assez gros échantillon, avec peu de loci très polymorphes.

129L’absence de tout déséquilibre de liaison est confirmée par les nouvelles analyses (p-value > 0.1).

130Il y a du stuttering très significatif à tous les loci, sauf IR08. La différence avec les analyses de la 1^re édition provient du fait que j’ai adopté la stratégie décrite dans mon article de 2019 (De Meeûs et al., 2019) où on regarde s’il y a un déficit d’hétérozygotes entre allèles distants d’une répétition et, pour les loci imparfaits (tous sauf IR27), ou deux répétitions sur les sorties graphiques de MicroChecker, et avec 10 000 randomisations. Cependant, une tentative de correction effectuée, en regroupant en allèles synthétiques les allèles proches en taille (De Meeûs et al., 2019), n’a abouti à aucune amélioration. La significativité provient sans doute du rôle important des allèles nuls et du fait que la plupart des allèles se suivent avec un seul pas de différence. Le stuttering serait donc ici artefactuel.

131J’ai retesté la dominance d’allèles courts sur l’ensemble des sous-échantillons (cf. fichier Excel sur mon site : http://www.t-de-meeus.fr/Data/DataLivreInitiation/IxodesResults.xlsx), comme décrit dans Manangwa et al. (2019). Il s’agit d’étudier la corrélation taille d’allèle/F_IT et aussi, en cas de doute, la régression pondérée par p_i (1 – p_i) (où pi est la fréquence de l’allèle) taille d’allèle/F_IT, ou même avec F_IS pour confirmer. Le locus IR27 présente toujours une dominance d’allèles courts, marginalement non significative avec le coefficient de Spearman (p-value = 0,0532), mais significative avec la régression pondérée (p-value = 0,0173, R² = 0,5276).

132Les données manquantes (allèles nuls possiblement homozygotes) expliquent 71 % de la variation du F_IS d’un locus à l’autre et 96 % si j’exclue IR27 (p-value < 0,0417 dans les deux cas avec la corrélation de Spearman). Ce qui écarte définitivement le rôle du stuttering. L’intercept de la régression correspondante (0 donnée manquante et donc allèles nuls en fréquence très faible), sans IR27, correspond à un F_IS = 0,0951. Cette valeur pourrait être expliquée par des croisements frère-sœur à hauteur de 35 % s’ils expliquaient la totalité de ce déficit, mais duquel on ne peut exclure une interaction avec un effet Wahlund entre entités faiblement différentiées. Comme suggéré dans la 1^re édition de ce manuel, une variance importante de survie d’une fratrie à l’autre, combinée à une agrégation résiduelle des tiques issues d’une même ponte, pourrait expliquer la part des déficits d’hétérozygotes non expliquée par les allèles nuls.

133Pour être vraiment rigoureux, il conviendrait d’éliminer le locus IR27 des autres analyses. Pour ce locus, 50 % du déficit semble être expliqué par la dominance des allèles courts et le reste par les mêmes facteurs que pour les autres loci. Il n’existe pas de méthode analytique pour corriger cet effet qui a un impact sur l’estimation du F_ST ou autres descripteurs de subdivision des populations. Nous verrons aussi qu’il est indispensable de corriger l’effet des allèles nuls avec la méthode de Chapuis et Estoup (2007) pour les mesures et tests de structuration, et notamment l’isolement par la distance (Séré et al., 2017). Ceci n’est cependant pas possible pour les tests de biais de structuration sexe ou pathogène spécifique (voir plus loin), ou pour les analyses hiérarchiques à plus de trois niveaux.

134J’ai refait le test de biais de structuration sexe-spécifique sans les loci IR08 et IR27. Les résultats restent comparables à ceux de la 1^re édition de ce manuel.

135Pour le test d’isolement par la distance en Suisse, pour limiter le nombre de tests, j’ai choisi de commencer avec la pente de la régression de Rousset et son intervalle de confiance de bootstraps (5 000), quitte à faire le test de Mantel avec la distance de Cavalli-Sforza et Edwards en cas de résultat ambiguë, pour gagner en puissance (voir Séré et al., 2017). J’ai utilisé pour ce faire FreeNA (Chapuis et Estoup, 2007) pour analyser les données recodées selon les recommandations du logiciel avec la correction ENA. Un tutoriel pas à pas pour les analyses FreeNA est disponible dans la page « Enseignements » de mon site web. Cette notice est en anglais mais reste assez facile à comprendre, même si vous parlez anglais comme un cétartiodactyle ibérique. J’ai effectué l’analyse sur les femelles seules et les cinq loci (avec 5 000 bootstraps) ou sur toutes les données mais sans le locus IR08 (et donc pas de bootstraps). J’ai récupéré les matrices des estimateurs de F_ST entre toutes les paires de sous-échantillons du fichier de sortie FreeNA « *.pFST » avec la correction ENA (la 2^e matrice) et ses intervalles de confiance de bootstrap (les deux dernières matrices) et n’ai conservé que les paires contemporaines (95 ou 96) pour construire un tableau de paires de sous-échantillons contemporains. J’ai ensuite calculé la pente de la régression et son intervalle de confiance avec Excel avec les fonctions graphiques « nuage de points », « courbes de tendances », « autres options », et « afficher l’équation sur le graphique » pour la pente moyenne et chacun de ses intervalles de confiance (IC 95 %). Les valeurs obtenues chez les femelles sont alors b = 0,0027 compris dans IC 95 % = [0,0013 ; 0,0052]. Puisque cet intervalle de confiance ne comprend pas le 0, l’isolement par la distance est donc significatif.

136Pour les effectifs efficaces, j’ai utilisé la dernière mouture de NeEstimator (Do et al., 2014) avec la méthode des déséquilibres de liaison (Waples et Do, 2008) corrigée pour données manquantes (Peel et al., 2013) et celle des co-ascendances de Nomura (Nomura, 2008). J’ai aussi utilisé Estim (Vitalis et Couvet, 2001 ; Vitalis, 2002). J’ai calculé les moyennes pondérées sur l’ensemble des valeurs obtenues en Suisse, les minimales et les maximales. Dans la feuille Excel mentionnée ci-dessus, on peut voir que la densité efficace des tiques (pour une surface d’échantillons d’environ 0,2 km² comme pour la 1^re édition de ce manuel) varie entre D_{e_min} = 1 203 et D_{e_max} = 6 100 tiques adultes par km² et une moyenne de 3 149 tiques par km², ce qui est équivalent aux résultats vus dans la 1^re édition et donc toujours cohérent avec ce que l’on observe sur le terrain. Il est important de rajouter que si la distribution des tiques est plus disparate que ce que nous avons rencontré dans nos sites, les valeurs obtenues représentent une surestimation des densités réelles des tiques en Suisse. Le calcul de la distance de dispersion par génération (en ordre de grandeur) est différent de ce que j’indiquais dans la 1^re édition, car je n’avais pas compris que σ est la distance axiale entre parents et descendants adultes et est égale à deux fois la distance de dispersion (δ) et que σ² de Rousset n’est pas le carré de cette distance mais la moyenne de ses carrés (qui aurait donc dû être symbolisée par afin d’éviter la confusion). Je remercie à ce titre Olivier Hardy de m’avoir fait remarquer ces erreurs lors d’une révision de l’article de Séré et al. (2017). L’erreur commise ne change pas l’ordre de grandeur des valeurs obtenues. Il faut simplement multiplier les anciennes valeurs par deux et ne pas oublier que si la variance de dispersion est grande la racine carrée de la moyenne des carrés ne représente qu’une estimation très approximative de la valeur moyenne réelle. La distance de dispersion approximative devient donc :

137Ici nous obtenons en moyenne δ ≈ 194 m avec une amplitude comprise entre 100 et 451 m par génération. Ces valeurs sont susceptibles d’être revues à la hausse en cas de surestimation des densités, mais ne devraient pas, en tout état de cause, dépasser les 4 km par génération. Ces valeurs sont très proches de celles que nous avions obtenues avec les clusters BAPS dans la 1^re édition. Les conclusions restent donc largement les mêmes : l’environnement est très visqueux pour les tiques qui ne peuvent envoyer leurs gènes à plus de 500 m (et au grand maximum 4 km) tous les 2-3 ans. J’en conclue que les paires d’individus de chaque cluster BAPS que nous avions gardées étaient un sous-échantillon représentatif de l’échantillon complet, et le fait d’avoir sélectionné les individus présentant un génotype complet avait probablement permis de minimiser l’effet des allèles nuls sur l’estimation de la pente (quand même plus élevée) et la puissance du test. Cela illustre aussi la nécessité de corriger pour les effets des allèles nuls. J’ai aussi fait la régression à partir du jeu de données complet sans IR08 et des femelles sans IR27. Avec quatre loci (toutes les données sans IR08, femelles sans IR27), il n’a pas été possible de réaliser de bootstraps. J’ai donc testé la significativité de cette régression par un test de Mantel. Comme les matrices de distances ne sont pas carrées, car il y a les paires pour 1995 et celles pour 1996, mais pas celles inter-années, j’ai dû le faire avec l’option « Mantelize it » de Fstat. Comme ce test est bilatéral, et que la pente était positive, j’ai divisé la p-value par 2. Les pentes restent similaires (b = 0,0028 et b = 0,0026 respectivement). Le test est significatif (p-value = 0,0059 et p-value = 0,0208 respectivement). Pour les données sans IR08, les effectifs efficaces sont légèrement plus faibles mais ne sortent significativement pas de la gamme de valeurs observées avec les femelles et tous les loci.

138Le schéma selon lequel il y a une forte variance de survie entre pontes et une sorte de signature spatiale des fratries (tendance même légère à rester ensemble) permet d’expliquer, en grande partie au moins, la structure génétique locale. Cette même variance de survie permettrait d’expliquer aussi le biais de structuration sexe spécifique si les sœurs sont davantage affectées que les mâles (mais on ne voit pas très bien par quel truchement). La dispersion des stades immatures s’effectue essentiellement probablement par les oiseaux. Dans ce cas, et surtout pour les larves, il est probable que de nombreux membres d’une même fratrie soient transportés ensemble. Cela permettrait de mieux comprendre pourquoi l’indice d’assignement donne un résultat significatif pour le test du biais de structuration sexe spécifique et pas le F_ST. Ici, pour expliquer le biais de structuration en faveur des femelles observé, il est possible que les femelles dispersantes survivent en moyenne très peu à ces voyages alors que les mâles persisteraient beaucoup plus facilement dans leur nouvel environnement. Mais une fois encore, il est difficile de comprendre par quel mécanisme. Une spécificité différentielle vis-à-vis des hôtes en fonction du sexe pourrait aboutir à ce signal si les larves et nymphes femelles délaissent les oiseaux au profit des mammifères moins mobiles, en particulier les micromammifères. Pour tester cette hypothèse, il faudrait pouvoir sexer les larves et nymphes trouvées sur différents types d’hôtes.

Interactions avec les micropathogènes transmis

Introduction

139La tique I. ricinus transmet un très grand nombre de pathogènes à ses multiples hôtes, dont la borréliose de Lyme qui, dans les régions boréales, représente un poids économique et en santé publique important (Gubler, 1998). Les agents de la borréliose de Lyme appartiennent au complexe d’espèces Borrelia burgdorferi sl. Il existe actuellement 21-22 espèces (ou génoespèces) dans ce complexe d’espèces, dont trois sont les agents majeurs de la maladie de Lyme : B. burgdorferi ss, B. afzelii et B. garinii ; sept sont des agents mineurs de la maladie ; dix sont de pathogénicité inconnue ; et une souche attend d’être mieux caractérisée (Eisen, 2020). Ces différentes espèces ne sont d’ailleurs pas responsables de symptômes identiques et présentent des spécificités d’hôtes réservoirs différentes (De Meeûs et al., 2004b). En Europe de l’Ouest, B. burgdorferi est préférentiellement retrouvée chez l’écureuil roux, B. afzelii chez des campagnoles, des mulots et aussi l’écureuil roux, B. garinii plutôt chez des oiseaux et B. spielmanii uniquement chez le loir (Richter et al., 2006 ; Postic et al., 2007). Il semble que le statut de B. burgdorferi ss (ss) ne soit pas aussi simple. Son association avec l’écureuil pourrait ne pas être absolue, même en Suisse où cette association semble prononcée (Humair et Gern, 1998). D’après la littérature récente, cette espèce serait généraliste sur la totalité de son aire de répartition mais plus spécialisée localement (Lin et al., 2020), en particulier en Europe continentale où elle montre une préférence pour les micromammifères (Berret et Voordouw, 2015) et peut-être pour l’écureuil en Suisse (Humair et Gern, 1998). L’épidémiologie de ces pathogènes reste largement mal connue et les résultats obtenus précédemment par nos analyses suscitent un certain nombre de questions. S’il y a spécificité différente des tiques immatures, sachant que les borrélies sont spécifiques des hôtes, les tiques des deux sexes devraient présenter des prévalences différentes pour les différentes espèces de borrélies. En particulier, les femelles devraient porter davantage de borrélies d’hôtes peu mobiles (B. burgdorferi, B. afzelii) et les mâles celles d’hôtes plus mobiles (B. garinii, B. valaisiana). Ensuite, il est possible que l’infection par les borrélies puisse modifier le schéma de migration. Enfin, dans la mesure où un conflit/coopération pourrait exister au sein des tiques, existe-t-il une corrélation entre la présence des différentes espèces de borrélies au sein de tiques ?

Présentation des données

140Toutes les tiques échantillonnées en Suisse pour cette étude avaient été coupées en deux, et une moitié envoyée à l’Institut de zoologie de Neuchâtel pour détermination de présence de borrélies et détermination de l’espèce (sondes moléculaires). L’autre moitié a été gardée dans l’alcool et un grand nombre utilisé pour génotypage microsatellite. Les données sont contenues dans le fichier TotBrutBorIR.txt où toutes les informations nécessaires sont disponibles. La présence ou l’absence de chaque espèce de borrélie trouvée est notée par un 1 ou un 0 dans la colonne correspondante. Un grand nombre de borrélies n’ont pu être déterminées au niveau de l’espèce (colonne “Bbundet”) et seules trois espèces ont été trouvées : B. burgdorferi (Bbss), B. afzelii (Bba) et B. garinii (Bbg, trouvée trois fois).

Distribution des différentes borrélies dans les femelles et mâles d’I. ricinus : analyses de la 1^re édition

141Pour cette analyse, nous allons devoir effectuer une régression logistique pour chaque espèce de borrélie (Bbundet, Bbss, et Bba). Bbg, trop rare sera laissée de côté. On va chercher à expliquer la présence de telle ou telle autre espèce de borrélie par le site, l’année et le sexe de la tique, ainsi que les interactions. Nous allons donc avoir besoin de R une fois de plus. Comme c’est le sexe que l’on souhaite tester ici, nous allons mettre ce facteur en premier (l’ordre compte dans les modèles de R). Après avoir lancé R et s’être positionné dans le répertoire approprié, on tape les commandes suivantes :

142afin de faire lire l’ensemble du jeu de données à R (NB le > est automatiquement inséré par R). On spécifie ensuite le modèle en tapant la commande (sur une ligne) :

143On remarque que l’interaction entre facteurs est codée avec un “:” et que la régression est logistique, car on spécifie bien qu’elle appartient à la famille binomiale avec un lien “logit” de la moyenne. Le lien logit signifie juste que la fonction qui relie la probabilité moyenne de la variable à expliquer (P_Bba probabilité de trouver une Bba) est du type log(P_Bba/(1 – P_Bba)) et la variance égale à P_Bba/(1 – P_Bba). Dans notre cas, la variance est en fait inférieure à cette valeur et il y a sous-dispersion, ce dont nous discuterons plus loin.

144Ensuite, il s’agit de tester le modèle par la commande :

145Le test est en effet un Chi2, car nous comparons des fréquences. Cette commande renvoie au résultat suivant :

146Nous constatons que seuls les termes “Site” et “Year” semblent importer et que le logiciel n’est apparemment pas très satisfait de la qualité des données. Pour simplifier ce modèle, une commande pratique est la commande “step” qui permet d’analyser la qualité de différents modèles plus simples en retirant et ajoutant des termes l’un après l’autre en commençant par les interactions d’ordre supérieur (celles faisant appel au plus grand nombre de facteurs). Ceci est évalué à l’aide d’un critère appelé AIC (Akaike Information Criterion) (Akaike, 1974) dont la valeur, qui doit être minimisée, est une mesure de la qualité d’ajustement du modèle statistique estimé par rapport aux données. Il ne s’agit pas d’un test, mais d’un outil d’aide à la sélection du modèle le plus simple permettant d’expliquer au mieux les données, le modèle doté du plus petit AIC étant le meilleur (cf. réponse 12 pour plus de précisions). En tapant donc la commande :

147nous obtenons les résultats pour une série de différents modèles de plus en plus simples où les différents termes sont retirés un à un en commençant par l’interaction la plus complexe (Sex:Site:Year), qui est éliminée, l’AIC obtenu (338,69) s’avérant inférieur à celui du modèle complet (343,81), puis les interactions plus simples (Sex:Site et Sex:Year), jusqu’à ce que le retrait des facteurs conduisent à une augmentation de l’AIC par rapport au précédent. Ci-dessous sont présentés le début et la fin du processus :

148La dernière ligne présentée ci-dessus donne le meilleur modèle. Suivent des informations sur les coefficients associés aux différents facteurs que nous n’allons pas utiliser, ainsi que des messages d’alertes sur la mauvaise qualité des données (on ne fait pas de miracles). Il s’agit maintenant d’analyser en détail ce meilleur modèle avec la série d’instructions (pour gagner du temps on peut copier le modèle ci-dessus et le coller après avoir tapé "Bba2<-") :

149qui renvoie au résultat :

150La conclusion est donc qu’en ce qui concerne Bba, seuls le site et l’année importent. Ils expliquent respectivement 100 × 35,72/358,68 = 10 % et 100 × 8,86/358,68 = 2 % de la déviance totale. En procédant d’une manière identique pour Bbg, nous observons qu’aucune des variables n’explique les données alors que pour Bbss, en plus du site qui explique 28 % de la déviance totale (P-value < 0,001), le sexe des tique explique 3 % de la déviance (P-value = 0,007). Enfin, pour Bbundet le site seul explique 15 % de la déviance totale (P-value < 0,001).

151Comme je l’ai déjà signalé plus haut, la dispersion des résidus ne suit probablement pas une loi binomiale et la variance est probablement différente de P/(1 – P). Pour vérifier cela, il faut calculer le paramètre φ = Var(µ) × (1 – µ)/µ qui est ici inférieur à 1 (sous-dispersion) en particulier pour Bbss. On peut le calculer facilement avec la fonction "quasibinomial" (voir réponse 13). Comme seul Bbss a donné quelque chose de significatif pour le sexe des tiques, nous allons vérifier cela sur cette bactérie. Sous R, après avoir chargé le fichier de données si ce n’est déjà fait, nous allons taper les instructions suivantes :

152ce qui renvoie au résultat suivant (je ne garde que ce qui est le plus utile) :

153Nous voyons donc que le le coefficient de dispersion est petit (0,52), il y a donc bien sous-dispersion (pour Bba Φ = 0,93, il n’y a pratiquement pas de sous-dispersion pour cette espèce-ci). Nous voyons également que le sexe des tiques est important (significatif) avec une estimation négative pour les mâles (les données partielles, corrigées des autres effets, sont centrées sur 0). Ceci est vérifiable en tapant la commande anova(Bbss, test="F") (les modèles quasi se testent avec un F), ce qui donne :

154Sachant que le comportement des modèles quasi en régression logistique peut s’avérer étrange quand l’événement étudié (présence de Bbss) est rare, ce qui est notre cas, on est en droit de chercher à renforcer ce résultat. En fin de compte, nous cherchons juste à vérifier si nous n’avons pas plus de Bbss chez les tiques femelles que chez les mâles, puisque ces borrélies sont spécifiques de petits rongeurs peu dispersants, supposés être davantage parasités par les larves et nymphes femelles que mâles, quel que soit le site ou l’année. On peut donc calculer parmi les tiques infectées par Bbss, la proportion de tiques femelles et mâles et comparer cette proportion à ½ par un test binomial. Sur 34 tiques infestées par Bbss, 26 étaient femelles, ce qui conduit à la P-value du test binomial (sous R, binom.test(26, 34, p=0.5, alternative="greater")) P_{bino_26/34,0.5} = 0,0015, ce qui est équivalent aux résultats précédents. Vous vous demandez alors pourquoi vous ai-je cassé les pieds avec toutes ces régressions, alors qu’il était si simple de commencer par le test binomial ? La réponse est simple. D’abord, il n’est pas inutile d’apprendre à taquiner les régressions linéaires généralisées qui servent très souvent et, ensuite, dans une publication, une régression logistique en « quasi-likelihood » va avoir beaucoup plus de classe (en apparence) qu’un petit test binomial et impressionner beaucoup plus facilement ces referees désobligeants qui empoisonnent si souvent nos soumissions d’articles.

155Donc Bbss, borrélie d’écureuil en Suisse, est plus fréquente chez les tiques adultes femelles que mâles, suggérant ainsi une prédisposition de ces femelles à se nourrir sur cet hôte quand elles sont aux stades larvaire et/ou nymphal.

Analyses correctes en modèles généralisés pour cette réédition

156Grâce à Renaud Lancelot je sais maintenant comment analyser les résultats d’un glm sous R sans que l’ordre d’entrée des variables importe. J’ai aussi évité d’utiliser la fonction quasibinomiale, car finalement, je ne suis pas certain que cela corrige les problèmes de dispersion des résidus. Les lignes de commande pour ce faire sont les suivantes :

157ou bien éventuellement pour gagner du temps, nous pouvons commencer par une procédure de sélection du modèle minimum avec la commande « step »:

158on peut maintenant analyser le modèle minimum.

159Seul le site importe donc avec ces analyses pour Bbss. Pour les autres espèces, le facteur sexe ne reste pas dans le modèle minimal.

160En fait, la distribution des résidus et le fait que le test soit bilatéral pourraient expliquer ces résultats négatifs et, effectivement, le test binomial unilatéral de la 1^re édition est sans doute ce qui convient le mieux. J’ai aussi effectué un test exact de Fisher avec Rcmdr (test bilatéral donc moins puissant) et les conclusions restent les mêmes : les femelles sont plus souvent trouvées infectées que les mâles.

161Ces résultats sont présentés dans la feuille de calcul “IxodesResults.xlsx”.

Co-occurrence des différentes espèces de borrélies

Analyses de la 1^re édition

162Les différentes espèces de borrélies peuvent se retrouver en compétition, car elles partagent la même espèce de vecteur. Elles sont donc potentiellement en conflit et on pourrait s’attendre à un évitement. Au contraire, il pourrait y avoir association positive si les intérêts convergent ou si l’une des deux espèces immunodéprime ses hôtes et favorise ainsi l’entrée d’autres pathogènes. Il est donc intéressant de tester si ces borrélies se rencontrent au hasard ou non. La problématique est identique à une recherche d’association statistique entre deux états (infecté/non infecté) de deux caractères (espèce x, espèce y). On peut donc simplement appliquer la même procédure que pour un test de déséquilibre de liaison. Il suffit donc de coder la présence de chaque borrélie comme un locus et l’absence par 11 et la présence par 22. Il y a donc quatre loci (Bba, Bbg, Bbss, Bbundet) avec chacun deux allèles (1 ou 2), toujours homozygotes (ou haploïdes). Pour ce faire, il suffit d’ouvrir le fichier “TotBrutBorIR.txt” et d’y remplacer, dans l’ordre, tous les 1 en 22 et tous les 0 en 11 et de fusionner les colonnes Site year sex pour obtenir quelque chose de la forme (fig. 43).

Figure 43

163Enregistrons ce fichier en le nommant “TotBrutBorIRCoOccur.txt” et importons-le dans Genetix afin de le convertir au format Fstat. Cliquez sur Fichier, Importer. Choisissez l’option fichier texte et double-cliquez sur “TotBrutBorIRCoOccur.txt”. Choisissez les options séparateur tabulation, un chiffre par allèle, décochez la case de l’identifiant des individus et cliquez sur OK.

164Il faut ensuite cliquer sur le menu Link. Dis et choisir Black & Kafsur comme sur la figure 44, ce qui aura pour effet de lancer une fenêtre de choix que vous devrez rendre comme dans la figure 45.

Figure 44

Figure 45

165Cliquez ensuite sur OK et les résultats s’affichent dans TotBrutBorIRCoOccur.lkd.

166Cliquez ensuite sur Outils, Conversion et FSTAT et nommez le fichier “TotBrutBorIRCoOccur.dat”. Genetix construira donc un fichier où seront considérées comme appartenant à des populations différentes les tiques de sites, d’années et de sexes différents. Ouvrons ce fichier sous Fstat et sélectionnons les mêmes options qu’en figure 46.

Figure 46

167Constatez que nous ne gardons que les fréquences alléliques (cela pourrait servir) et ce qui nous intéresse, le test de déséquilibre de liaison. On choisit dans un premier temps le niveau 5/100 pour aller plus vite. Cliquez sur “Run” et ensuite ouvrez le fichier “TotBrutBorIRCoOccur.out”. Vous constatez que seulement 2 640 permutations ont été effectuées. Recommencez donc avec le niveau 1/100 pour le menu “Nominal level for multiple tests”. Le résultat peut être synthétisé dans le tableau 18. On y voit clairement une association positive entre Bbss, Bba et Bbg, même si les associations avec Bbg sont marginalement significatives, on peut considérer que le signal existe eu égard à la grande rareté de Bbg (puissance très faible du test). Il est intéressant de noter pour information que Bbundet, vraisemblablement composée d’une mixture de Bbg (très largement sous-représentée ici) et Bbv (B. valaisiana curieusement absente de l’échantillon) donnent des valeurs essentiellement négatives pour R(IJ), ce qui fait regretter plus encore que les déterminations de l’époque aient connu autant de problèmes. Il n’en reste pas moins qu’une forte corrélation positive lie Bbss, Bbg et Bba, qui est confirmée si on teste la co-occurrence des trois espèces dans la même tique rencontrée une fois dans l’échantillon des 73 tiques femelles de Neuchâtel en 1996, et pas à Bern comme annoncé dans la partie résultat de l’article de De Meeûs et al. (2004b) (on ne relit jamais assez ses épreuves). Il y a N = 73 observations, une fréquence observée de 4/73, 6/73 et 1/73 pour Bbss, Bbg et Bba respectivement, donc une fréquence attendue de p = (4 × 6 × 1)/(73)³ pour l’événement de co-occurrence des trois borrélies dans la même tique, événement observé avec la fréquence k = 1. Cette fréquence observée peut être comparée à l’attendue par un test binomial. Sous R, tapez “binom.test(1, 73, p=0.00006169, alternative="two.sided")”, ce qui donne une P-value = 0,0045 très significative. Cette P-value est en fait égale à la probabilité de l’événement lui-même puisqu’il n’y en a pas de plus rare possible. Elle est donc égale à la probabilité (dans une loi binomiale) de tirer une seule fois Bbss-Bbg-Bba dans 73 tirages et où la probabilité de tirer Bbss-Bbg-Bba une fois (un tirage aléatoire) est de 0,00006169, soit (cf. n’importe quel manuel de statistiques) :

Tableau 18 – Valeurs (moyennes pondérées par les tailles de sous-échantillons) de corrélation entre la présence de chacune des deux bactéries considérées (ligne 2-2 dans la sortie de Genetix, colonne R(IJ)) et P-value (test G multi-échantillon de Fstat) correspondante. Une valeur de R(IJ) positive indique une association positive, alors qu’une valeur négative indique une répulsion.

168Cette corrélation est donc très forte. Elle peut être due au fait que les tiques infectées correspondent à des individus sensibles et que les autres individus sont résistants. Cette corrélation peut également provenir du fait qu’être infecté par une des trois borrélies tend à favoriser l’infection par les deux autres (par immunosuppression, par exemple). Ceci peut être testé en ne regardant que les tiques infectées. La corrélation existe-t-elle toujours ? Ici j’ai dû effacer les analyses initiales qui donnaient en fait des résultats biaisés. Nous allons donc directement consulter les résultats des analyses effectuées pour la 2^e édition de ce livre.

Analyses effectuées pour la réédition de ce manuel sur les occurrences des différentes espèces de borrélies dans leur environnement

Co-occurrences des différentes espèces de borrélies

169Ici, il convient d’utiliser une autre technique que celle décrite dans la 1^re version et qui donne en fait des résultats biaisés quand on s’adresse aux tiques infectées (saines exclues). J’ai utilisé les abréviations ss, a, g et undet pour, sensu stricto, afzelii, garinii et non-déterminés, respectivement. J’ai effectué des tests binomiaux sous R pour comparer la fréquence des occurrences observées à celles attendues (produits des fréquences de chaque acteur), sur l’ensemble des données, et sur les tiques infectées seules (tiques saines éliminées). Avec Nss, Na et Nundet les nombres de tiques infectées par les différentes espèces de borrélies NObs le nombre de co-occurrences observées (11 pour ss et a), NTot le nombre total d’occurrences (ici 858 en comptant toutes les données), Pss, Pa, et Pundet les prévalences des différentes borrélies (Pss = Nss/NTot), et PropExp la proportion attendue de l’occurrence (PropExp = Pss × Pa = 0,0396 × 0,0536 = 0,002 pour ss avec a), la commande sous R était donc du type : binom.test (NObs, NTot, PropExp). Les résultats sont compilés dans le tableau 19.

170Sur l’ensemble des données, les résultats sont quasiment les mêmes que dans la 1^re édition, avec une tendance confirmée d’association positive entre espèces de borrélies.

171En revanche, pour les données sans les tiques saines (indemnes de Borrélies), les résultats diffèrent quelque peu et sont d’ailleurs rendus plus faciles à interpréter biologiquement. En effet, il n’y a plus exclusion significative (valeur négative pour DR dans le tableau 19) entre ss et aa qui partagent en effet le même type d’hôtes vertébrés (micromammifères), dans la zone d’étude. Par contre, pour les borrélies non déterminées, que nous pouvons largement soupçonner d’appartenir aux génoespèces d’oiseaux (B. garinii et B. valaisiana), ces dernières semblent éviter de se retrouver trop souvent avec des borrélies de micromammifères (ss et a), avec DR négatifs et très significatifs (tabl. 19).

172Les détails des analyses sont présentés dans la feuille de calcul « IxodesResults.xlsx ».

173Ces résultats cadrent bien avec une vulnérabilité plus grande de certaines tiques vis-à-vis des borrélies dans leur ensemble. Ceci crée une co-occurrence positive artificielle. Il n’est pas vraiment dans l’intérêt des différentes espèces de borrélies qui sont spécifiques d’hôtes différents de se retrouver dans une même tique qui ne va se nourrir que sur un seul type d’hôte à la fois et ne pouvoir transmettre qu’une seule espèce. Seules ss et a, spécifiques de micromammifères, peuvent être indifférentes à cela, si tant est que la tique où ils cohabitent se nourrisse sur un écureuil, au cas où la spécificité des ss suisses vis-à-vis de cet hôte serait confirmée (ce qui n’est en fait pas vraiment le cas). Il est donc cohérent, parmi les tiques réceptives à ces borrélies, qu’il y ait évitement entre borrélies d’oiseaux et borrélies de micromammifères. Comme discuté dans la 1^re édition, cela résulte soit d’une adéquation de nature gène pour gène entre génoespèces de borrélies et génotype de la tique, soit il y a tendance à l’exclusion d’un type de génoespèce par l’autre (par destruction directe ou médiée par le système de défense de la tique), soit par manipulation de la tique par la bactérie qui la parasite (favorisation), soit par simple contrainte spatiale : à l’échelle du site (un repas sur un type d’hôte donné au stade larve augmente la probabilité que la nymphe se nourrisse sur le même type d’hôte) ; ou à l’échelle globale (chaque site est colonisé plutôt par un seul type de borrélie). Comme par définition les borrélies d’oiseaux se trouvent chez les oiseaux et celles de micromammifères sur des micromammifères, pour qu’une tique se trouve infectée par les deux types de borrélies, il faut que cette dernière ait pris deux repas différents, sur deux types d’hôtes différents. Comme il n’y a pratiquement pas de transmission transovarienne pour ces borrélies (Norte et al., 2020), la présence dans la même tique adulte de spirochètes d’oiseaux et de micromammifères signifie un repas sanguin chez l’un pour la larve et chez l’autre pour la nymphe ou inversement.

Distribution spatiale

174Nous nous focaliserons sur les trois génoespèces Bbss, Ba et Bbundet (pas assez de Bg) et pour l’année 1996 (pas assez de données en 1995). Les résultats discutés sont illustrés par la feuille de calcul “IxodesResults.xlsx”.

175Ces trois espèces présentent globalement une structure agrégée sur les sites, comme l’indiquent des rapports variance sur moyenne d’occurrences sur sites très supérieurs à 1 (attendu sous l’hypothèse de distribution aléatoire de Poisson). Cela signifie que les borrélies sont surtout présentes dans de rares sites et plutôt absentes partout ailleurs. Elles ne sont pas présentes dans les mêmes sites, comme le montre un test du Chi² (p-value = 0,0002), et la p-value = 0,0165 combinée des tests exacts de Fisher exécutés par site (huit tests), que j’ai effectué en plus afin de vérifier la validité du Chi² (six classes attendues se montrent inférieures à 1 parmi les 48). Donc chaque espèce de borrélie occupe des espaces qui lui sont plutôt propres. Cela corrobore l’hypothèse de contrainte spatiale pour expliquer l’exclusion des espèces de borrélies entre elles.

176Enfin, des tests exacts de Fisher n’ont pas permis de voir de différences de prévalences globales entre les trois espèces (p-value > 0,1) ni de sa variance (test de Siegel-Tuckey, Siegel et Castellane, 1988 ; p-value > 0.058). Tous ces tests ont été effectués avec le package Rcmdr de R.

177Personnellement, je ne comprends pas très bien pourquoi certains sites seraient plus favorables à certaines borrélies et moins à d’autres, étant donné qu’à priori les mêmes cohortes d’hôtes vertébrés y cohabitent. Les borrélies d’oiseaux, en particulier, devraient pouvoir se trouver n’importe où. Peut-être faut-il y voir une contrainte historique et spatiale qui fait que certaines colonies de tiques se sont retrouvées inféodées à des oiseaux et à leurs microsites de visites qu’ils ne partagent que peu avec les micromammifères, alors que d’autres colonies de tiques se retrouvent établies dans des microsites à micromammifères.

Occurrence des différentes espèces de borrélies et génétique des tiques : analyses de la 1^re édition

178Dans cette partie, nous rechercherons s’il existe une relation entre la génétique des tiques et leur probabilité d’infection par chacun des quatre types de borrélies. On peut répondre à cette question de trois manières. Soit en testant la différenciation génétique entre tiques infectées et non infectées dans chaque sous-échantillon, soit en testant la différenciation, dans chaque sous-échantillon, entre tiques infectées par des borrélies différentes, enfin en procédant à un test de biais de structuration, comme nous l’avons fait pour le sexe des tiques, mais avec le statut infecté/non infecté à la place.

Différenciation entre tiques infectées et non infectées

179Il faut construire un fichier par espèce de bactérie Bbss, Bba et Bbundet (il n’y a pas assez de Bbg). On doit changer de nom de population pour chaque site, année et sexe. La figure 47 donne un exemple de fichier pour Bbss.

Figure 47

180Il suffit ensuite de convertir ce fichier au format Fstat (en passant par Genetix, par exemple) et de procéder sous Fstat au calcul des F_ST par paire de sous-échantillons et au test de différenciation par paire, comme indiqué dans la figure 48. Vous constatez que j’ai coché la case 1/1000 pour le nominal level afin d’obtenir au moins 10 000 permutations et donc d’obtenir des P-values assez précises. Le fichier de données s’appelle “ForPairedBbss.dat” et les fichiers de sortie qui nous intéressent sont “ForPairedBbss.fst” pour récupérer les valeurs de F_ST par paire qui nous intéressent et “ForPairedBbss-pp.pvl” où nous allons récupérer les P-values correspondantes. Attention, dans ces fichiers, seules les comparaisons entre tiques infectées et non infectées du même sexe, de la même année et du même site nous intéressent. Le résultat pour Bbss est présenté dans le tableau 20.

Figure 48

181Vous remarquerez que la combinaison des cinq tests disponibles a été effectuée à l’aide de la procédure binomiale généralisée de Teriokhin et al. (2007) effectuée à l’aide du logiciel MultiTest (De Meeûs et al., 2009). En effet, à partir de quatre tests, je préfère utiliser cette procédure plutôt que le test Z de Stouffer (Whitlock, 2005). Pour k < 4, nous pouvons aussi utiliser la procédure binomiale mais avec k’ = k (De Meeûs, 2004). Pour effectuer le Z de Stouffer, chaque P-value individuelle est transformée en son équivalent de la distribution Z centrée sur 0 et d’écart-type 1. Sous Excel, on tape =SI(B2="NA";"";SI(B2>0.9999;LOI.NORMALE.INVERSE(0.9999;0;1);LOI.NORMALE.INVERSE(B2;0;1))). B2 correspond aux coordonnées de la case du tableau Excel où la P-value à transformer se trouve. Cette commande renvoie une absence de résultat quand “NA” est rencontré et tient compte du fait qu’une P-value de 1 n’est pas transformable et la P-value = 0,9999 est choisie comme limite supérieure. Enfin, l’équivalent de la P-value en Z centrée réduite de moyenne 0 et d’écart-type 1 est calculé. Les valeurs Z_i obtenues sont ensuite combinées dans la formule (Whitlock, 2005) :

182 , où k est le nombre de tests (= 0,3266 ici).

183La P-value globale s’obtient ensuite par un retour à la loi normale, soit sous Excel : =LOI.NORMALE.STANDARD(Z_s)(=0,628 ici). Vous trouverez un argumentaire plus détaillé dans De Meeûs et al. (2009) pour les situations où la procédure binomiale généralisée ou le test Z doivent ou peuvent être utilisés.

184Si on procède de la même façon pour Bba et Bbundet, le même type de résultat émerge, même quand on ne distingue pas le sexe des tiques (échantillons plus grands) puisque pour ces deux catégories de borrélies, nous avons vu que le sexe des tiques n’importait pas. Ce résultat est rassurant car, étant donné que les marqueurs sont non codants (donc neutres) et indépendants, il eut été difficile d’interpréter une différenciation entre tiques infectées et non infectées, à moins d’évoquer l’existence d’espèces cryptiques de tiques et une spécificité des borrélies.

Différenciation entre tiques infectées par différentes borrélies

185Ici, il faut ne garder que les tiques infectées et définir comme sous-population les tiques du même sexe, échantillonnées la même année, dans le même site et ayant le même statut infectieux. Notons qu’une tique infectée par Bba et Bbss ne fera pas partie de la même sous-population qu’une tique infectée par Bbss seule. On met ensuite le fichier au format Fstat et on lance la procédure de F_ST par paire. Ce faisant, vous constaterez que la plupart des tests sont infaisables, c’est normal. Les résultats sont compilés dans le tableau 21. En toute rigueur les tests, qui ne sont pas tous indépendants, devraient subir la correction de Bonferroni. Cependant, étant donné la faiblesse des échantillons (manque total de puissance), nous nous abstiendrons de le faire. Le seul F_ST positif est obtenu entre Bba et Bbundet, mais il n’est pas significativement plus grand que 0. Eu égard à la faiblesse des tailles de sous-populations ici, nous décidons que rien ne permet d’affirmer l’existence d’une différence génétique entre tiques infectées par différentes bactéries et rien ne permet de l’exclure formellement au moins pour ce qui concerne le couple Bba/Bbundet. S’il existe des races d’hôtes chez I. ricinus, ce n’est pas avec ces données qu’on peut le montrer.

Biais de structuration spécifique associé au pathogène

186Ici, il faut reprendre les données pour chaque espèce de bactérie et créer un fichier de type Genepop comme ce qui a été fait en p. 153-156, sauf qu’ici les tiques sont distinguées en fonction de leur statut infectieux et non par leur sexe, tel que dans la figure 49. Notons que nous ne traitons que les sites prélevés en Suisse et où au moins une tique infectée est trouvée. Parce qu’il y a un biais de structuration sexe-spécifique, ainsi que des différences d’infection, les femelles et les mâles sont analysés séparément. Cependant, parce que la taille des échantillons est très faible (peu de borrélies trouvées et identifiées), nous combinerons le tout dans un seul fichier (gain de puissance). On prendra soin de distinguer les tiques d’années et de sexe différents comme appartenant à des populations différentes (séparées par un “pop” dans le fichier).

Figure 49

187Quand le fichier est constitué, il faut ensuite lancer Fstat et cliquer sur le menu “Biased dispersal”. La fiche correspondante apparaît alors. Il faut ensuite charger le fichier à analyser en cliquant le menu “File” et “Open” et cocher les cases comme en figure 50 puis sur le bouton “Go!”. Pour une raison que j’ignore, il faut cocher tous les paramètres si on souhaite obtenir le résultat du test sur H_s, en particulier F_IS et H_o qui ne sont guères utiles ici, car nous avons codé les mâles homozygotes pour IR08.

Figure 50

188Le résultat est contenu dans un fichier de type nomdufichier.res (un fichier par espèce de borrélie). Le résultat principal concerne le test du F_ST (et aussi la relatedness, ce qui est normal si on regarde sa définition dans la documentation de Fstat) et est présenté dans le tableau 22.

189Il y a donc bien un biais de structuration dû à l’infection par Bba. Plusieurs hypothèses peuvent expliquer ce résultat. La première hypothèse implique que certaines tiques, plus sensibles à l’infection par Bba sont aussi pléiotropiquement moins mobiles. Les marqueurs utilisés étant des microsatellites non codants, cela impliquerait un déterminisme génomique peu vraisemblable. Par ailleurs, l’absence totale de différenciation entre tiques infectées et non infectées (montrée en p. 195-198) discrédite cette interprétation.

190La deuxième hypothèse implique l’existence d’au moins deux espèces cryptiques dont l’une, moins mobile que la seconde, serait plus sensible à l’infection par Bba. Notons que nous n’avons noté aucun déséquilibre de liaison (attendu en pareil cas). Par ailleurs, si on calcule avec Fstat le F_IS des tiques en séparant celles infectées par Bba de celles qui ne le sont pas, on ne retrouve aucune diminution du F_IS (~ 0,45 pour les infectées et ~ 0,44 pour les non infectées par Bba). Il n’existe pas de différenciation significative entre tiques infectées et non infectées. Cette interprétation n’est donc aucunement soutenue.

191La troisième interprétation possible impliquerait l’existence d’une adaptation locale des borrélies qui infecteraient plus facilement les tiques locales (résidentes) que les immigrantes. Deux arguments vont à l’encontre de cette hypothèse. La première est que les tiques mâles et femelles qui en principe n’ont pas la même dispersion (les femelles dispersent en principe peu ou pas, cf. p. 153-159) ne sont pas infectées différemment par Bba (p. 183). Par ailleurs, c’est le partenaire le plus mobile des deux qui doit en théorie être le mieux adapté localement (Gandon et al., 1996 ; Gandon, 2002). Or ici, les tiques sont modestement structurées alors que l’on pense que les borrélies le sont beaucoup plus (Qiu et al., 1997). C’est donc l’hôte (la tique) qui devrait être adapté localement et non l’inverse.

192La quatrième hypothèse implique une survie plus faible des tiques migrantes quand ces dernières sont infectées par Bba. Comme les tiques femelles sont moins mobiles que les mâles, ce sont ces derniers qui devraient être les plus affectés par ce phénomène. Ceci est testable en refaisant l’analyse sur les tiques femelles et mâles séparément. Cela suppose une survie au stress moins bonne des larves et/ou nymphes infectées par Bba.

193La cinquième hypothèse est la plus séduisante. Elle implique une manipulation des larves et nymphes par la borrélie. Cette borrélie est spécifique de petits rongeurs. Il est donc plus intéressant pour elle d’être injectée dans un petit rongeur, peu dispersant, que dans un oiseau ou un grand mammifère, hôtes beaucoup plus mobiles. Les Bba capables de manipuler les tiques qu’elles infectent de sorte que ces dernières préfèrent se fixer sur un petit rongeur plutôt que sur d’autres hôtes seraient donc avantagées. Cette hypothèse est testable en laboratoire, mais cela n’a malheureusement jamais été fait. Cela implique aussi, comme pour l’hypothèse précédente, que les femelles, déjà très peu mobiles, seront moins affectées par le biais de structuration Bba-spécifique que les mâles.

Biais de structuration spécifique au pathogène et au sexe

194Nous allons utiliser la même procédure que précédemment, mais en divisant le fichier en deux : un fichier pour les tiques femelles et un autre pour les tiques mâles. Cette fois, les tests seront faits de manière unilatérale avec I (infectés) comme catégorie la plus philopatrique. Il y a deux raisons à cela. La première est que l’on connaît d’avance le sens du signal. La seconde raison est que les échantillons étant encore plus petits, nous aurons besoin d’encore plus de puissance dans le test. Nous ne nous occuperons que du test sur le F_ST. Les tests sont tous les deux significatifs avec P-value = 0,0497 pour les tiques femelles et P-value = 0,0123 pour les tiques mâles et une apparente très forte différence de signal entre les deux, comme indiqué dans la figure 51.

Figure 51

195Nous pouvons également constater la formidable variance chez les mâles infectés (très peu nombreux). Nous pouvons effectuer un test unilatéral de Wilcoxon pour données appariées comme en p. 149 quand nous avions comparé les F_IS des données brutes avec ceux des données clusterisées par BAPS. Ici, l’unité d’appariement reste le locus (donc cinq données), mais la statistique est la différence de F_ST entre tiques infectées et non infectées chez les femelles et les mâles. Le fichier à tester contiendra donc les différences des différences appariées : (F_STMI – F_STMU) – (F_STFI – F_STFU). Le test unilatéral (F_STMI – F_STMU > F_STFI – F_STFU) montre que la différence n’est pas significative, même si la P-value reste relativement faible (0,17). Ceci illustre les limites de notre jeu de données (beaucoup trop de données manquantes).

Occurrence des différentes espèces de borrélies et génétique des tiques : nouvelles analyses

Différenciation entre tiques infectées et non infectées

196Avec Fstat, pour les femelles seules et tous les loci, ou pour toutes les données sans le locus IR08, aucun test de différenciation n’a été significatif (p-values > 0,16).

197Par contre, quand nous utilisons FreeNA avec les femelles (pour avoir cinq loci et les résultats de 5 000 bootstraps), nous constatons une différenciation très significative entre tiques infectées par une espèce de borrélie et celles non touchées par cette bactérie. Le F_ST moyen est de F_ST = 0,0105, celui corrigé par FreeNA F_ST-FreeNA = 0,0699 avec un intervalle de confiance à 95 % (des moyennes) IC 95 % = [0,0001, 0,1515] qui ne contient pas le 0. Nous observons dix intervalles de confiance ne contenant pas le 0 (significatifs à 5 %) sur 18, ce qui est très supérieur aux 5 % attendus sous l’hypothèse nulle (p-value < 0,0001). Il ne semble pas que le fait d’être infectées par telle ou telle autre espèce de borrélie ait une influence significative sur le niveau de différenciation d’avec les tiques non infectées.

Différenciation entre tiques infectées par différentes borrélies

198Ici encore, avec Fstat (donc sans corriger pour les allèles nuls), que ce soit pour les femelles et cinq loci ou toutes les données sur quatre loci (IR08 exclu), la plupart des F_ST sont faibles, voire négatifs, et les tests de différenciation non significatifs.

199En utilisant la correction pour les allèles nuls sur les femelles et leur cinq loci, nous ne trouvons pas de différenciation significative entre tiques infectées par des espèces différentes. Il existe par contre une signature significative de subdivision entre tiques saines de toute borrélie et celles infectées par Bbss, Bbundet ou par au moins deux espèces. Par ailleurs, les échantillons de tiques infectées étaient très modestes (manque de puissance), ce qui peut expliquer l’absence de tests significatifs entre tiques infectées par des espèces de borrélies différentes. Il n’en reste pas moins qu’avec les moyennes calculées avec les F_ST corrigés par FreeNA, on obtient une valeur positive, quelles que soient les entités comparées, avec des IC 95 % significatifs dès que l’on oppose les tiques saines aux tiques infectées.

200Combiner les tiques infectées par n’importe quelle borrélie et les opposer aux tiques saines aboutit à des résultats ambigus avec F_ST-FReeNA = 0,0184 et un IC 95 % = [– 0,0047, 0,0475], mais deux sites sur six avec une différenciation significative (p-value = 0,0328, test binomial) chez les femelles. Cela est confirmé par des tests de comparaisons entre groupes de Fstat où le F_ST est en effet toujours supérieur à 0 mais n’est jamais significatif. Le fait que l’on trouve une moyenne inférieure quand on oppose les tiques infectées par n’importe quelles borrélies aux tiques saines, par rapport aux niveaux de subdivision observés entre tiques non infectées et tiques infectées par une espèce en particulier, suggère fortement que les génotypes des tiques infectées par une espèce de borrélie particulière diffèrent légèrement de celles infectées par une autre espèce. En effet, le fait d’opposer tiques infectées en général aux tiques saines aurait dû conduire à une valeur similaire à la moyenne obtenue avec les infections par chaque espèce de borrélie.

Différenciation génétique pathogène spécifique

Tests de comparaisons de niveaux de subdivision

201Ici, j’ai considéré uniquement les données de 1996. J’ai effectué deux types de tests : j’ai comparé les paramètres de structure génétique entre tiques infectées par telle ou telle autre espèce de borrélies et celles non infectées par cette espèce en particulier et j’ai ensuite opposé les tiques infectées par n’importe quelle borrélie aux tiques saines. Je ne retranscrirai ici que les résultats obtenus avec le F_ST corrigé pour les allèles nuls avec FreeNA.

Pour Bbss

202Pour les femelles et les cinq loci, nous observons que les tiques infectées par Bbss (F_ST-FreeNA = 0,0059, IC 95 % = [– 0,0157, 0,0253]) ne sont pas plus ou moins subdivisées que celles non infectées (F_ST-FreeNA = 0,0021, IC 95 % = [– 0,0039, 0,0078]).

203Pour toutes les données sans IR08, bien que nous ne disposions pas d’intervalle de confiance (pas de bootstrap avec seulement 4 loci), les valeurs observées chez les tiques infectées (F_ST-FreeNA = 0,0198, MiniMax = [– 0,0201, 0,0214]) chevauchent largement celles des tiques non infectées par cette espèce (F_ST-FreeNA = 0,0035, MiniMax = [– 0,0014, 0,0112]) (MiniMax = gamme des valeurs obtenues).

Pour Ba

204Ici, les femelles infectées (F_ST-FreeNA = 0,0848, IC 95 % = [0,0304, 0,1394]) montrent une subdivision significativement plus grande que 0 et plus importante que les tiques non infectées par cette espèce de spirochète (F_ST-FreeNA = 0,0036, IC 95 % = [– 0,0027, 0,0112]). Seul le locus IR08 montre un F_ST-FreeNA < 0.

205Pour toutes les données sans IR08, le chevauchement ne concerne que IR39 (0,0085) chez les infectées et IR32 (0,0087) pour les non infectées avec F_ST-FreeNA = 0,0481, MiniMax = [0,0085, 0,0891] et F_ST-FreeNA = 0,004, MiniMax = [0,0007, 0,0087] respectivement.

Pour Bbundet

206Les femelles infectées par ce que nous supposons (avec peu de doutes) être les borrélies d’oiseaux présentent la plus forte subdivision observée (F_ST-FreeNA = 0,1199, IC 95 % = [0,0755, 0,196]), très au-dessus de 0 et de la valeur observée pour les tiques non infectées par cette bactérie (F_ST-FreeNA = 0,002, IC 95 % = [– 0,0016, 0,0055]). Il n’y a ici aucun recouvrement entre aucun locus (voir la figure dans la feuille Excel « IxodesResults.xlsx », feuille « ComparDifInfUninf » que l’on trouve sur mon site : http://www.t-de-meeus.fr/Data/DataLivreInitiation/IxodesResults.xlsx).

207Pour les données sans IR08, aucun chevauchement n’est observé entre tiques infectées (F_ST-FreeNA = 0,1203, MiniMax = [0,0489, 0,1994]) et les tiques indemnes de cette borrélie (F_ST-FreeNA = 0,0031, MiniMax = [0, 0,0055]).

Toutes borrélies confondues

208Le signal reste significatif, même si la différence est moins importante entre femelles (5 loci) avec un F_ST-FreeNA = 0,0218 et un IC 95 % = [0,0139, 0,0336] pour les infectées et F_ST-FreeNA = 0,0029 et un IC 95 % = [– 0,0004, 0,0054] pour les tiques indemnes de toute borrélie. Le niveau de différenciation est significativement inférieur à celui observé pour les tiques infectées par Ba et Bbundet, ce qui suggère que les unes et les autres ne sont pas complètement équivalentes (races ?).

209Pour les données totales sur quatre loci, j’ai préféré séparer les mâles des femelles. Il n’y a pas de chevauchement chez les femelles avec un F_ST-FreeNA = 0,0262 et un MiniMax = [0,015, 0,0454], alors que les valeurs obtenues pour les mâles infectés, avec F_ST-FreeNA = 0,0139 et un MiniMax = [– 0,0111, 0,0466] embrassent entièrement la gamme de valeurs des mâles sains (F_ST-FreeNA = 0,0047 et un IC 95 % = [0,0017, 0,0073].

Tests de biais de dispersion pathogène spécifique

210Ici, que ce soit pour les femelles avec cinq loci ou les données sans IR08, aucun test n’est paru significatif (p-value > 0,2). Je n’ai pas séparé mâles et femelles (données avec 4 loci) car les données devenaient trop fragmentaires.

Conclusion sur le statut infectieux et la génétique des tiques

211Il semble bien y avoir une différence génétique, même légère, entre tiques saines et tiques infectées par une espèce particulière de borrélie, et peut être aussi entre tiques infectées par des borrélies différentes. S’il existe des tiques préférant se nourrir sur des oiseaux et d’autres sur des micromammifères, comme cela est suggéré dans l’article de Kempf et al. (2011), cela risque de se refléter sur leur statut infectieux. Comme nous pouvons le vérifier sur le graphique de la feuille Excel « DifSpDif » du fichier « IxodesResults.xlsx » sur mon site : http://www.t-de-meeus.fr/Data/DataLivreInitiation/IxodesResults.xlsx, les niveaux de subdivision observés sont en effet comparables à la contribution moyenne de l’espèce hôte à la mesure de subdivision des tiques qui avait été observée dans l’article de Kempf et al. (2011). Avec un déficit local d’hétérozygotes F_IS ≈ 0,1 (voir plus haut), si on accepte une contribution de la subdivision par l’hôte de 0,02 environ, il reste une valeur de F_IS ≈ 0,08 expliquée ni par les allèles nuls ni par les hôtes. Des croisements systématiques entre apparentés b = 0,3 pourraient expliquer cette valeur, soit par attractivité entre races d’hôtes, comme cela a été montré dans certaines populations (Kempf et al., 2009) soit par effets d’agrégations de fratries sur un même hôte des nymphes et surtout des larves et une variance importante de survie de ces fratries, qui pourraient aboutir à des différenciations génétiques associées au portage de certains types de borrélies. Une combinatoire des deux phénomènes représente peut-être l’interprétation la plus satisfaisante de ces résultats.

212Une autre conclusion est que la correction de l’effet des allèles nuls rend les résultats nettement plus clairs. Il existe donc bien un biais de structuration pathogène spécifique pour les tiques infectées par B. afzelii (comme montré auparavant) de micromammifères, mais aussi pour les indéterminées (potentiellement B.  garinii et B. valaisiana d’oiseaux). L’attraction préférentielle du type d’hôte ne semble plus trop convenir pour interpréter ces résultats dans la mesure où les oiseaux figurent parmi les hôtes susceptibles de disperser le mieux les tiques qui se trouvent dessus. Il semble donc qu’une survie plus faible des tiques dispersantes quand elles sont infectées reste l’hypothèse la plus vraisemblable pour expliquer nos observations.

Conclusions de la 1^re édition sur les borrélies et I. ricinus en Suisse

213Au cours de nos analyses, nous avons constaté que Bbss, borrélie d’écureuil, était plus souvent retrouvée chez les tiques mâles que femelles, ce qui est attendu si, comme le suggérait le biais de dispersion sexe-spécifique détecté chez ces tiques, les larves et nymphes femelles préfèrent se nourrir sur des rongeurs (peu dispersants). Rien de tel n’a pu être trouvé pour Bba pour laquelle ceci était attendu également, peut-être parce qu’une certaine quantité de tiques infectées par cette borrélie fait partie du stock Bbundet. Quant à Bbg, trop rarement détectée, d’autres études seront requises afin de déterminer si, comme attendu, elle est plus souvent retrouvée chez les tiques mâles.

214Certaines tiques sont plus sensibles ou plus exposées à l’infection par les borrélies en général, comme l’attestent les fortes corrélations positives observées sur les co-occurrences des trois espèces Bbss, Bba et Bbg. En se concentrant sur ces tiques sensibles (infectées par au moins une borrélie), il y a un évitement manifeste. Les corrélations deviennent toutes négatives, exception faite de l’association Bbss × Bbg, pour qui le faible nombre de Bbg détectées rend les choses difficiles à interpréter, et très significatives pour les couples Bba × Bbundet et Bbss × Bbundet. Cette dernière observation peut laisser spéculer que ces borrélies indéterminées soient majoritairement des borrélies d’oiseaux (Bbg et Bbv) très déficitaires dans notre jeu de données. Dans ce cas, nous pourrions proposer que les larves et nymphes sensibles se subdivisent en tiques ne se nourrissant que sur une gamme limitée d’hôtes réservoirs de borrélies spécifiques. Tout dépend de l’identité spécifique de ces Bbundet. Les données ne permettent pas d’exclure l’existence d’une telle spécificité en races d’hôtes. La manipulation de la spécificité des tiques par les borrélies ne peut pas non plus être exclue. C’est aussi cette manipulation qui expliquerait le biais de structuration des tiques infectées par Bba. D’une manière générale, on ne peut que regretter le nombre de données manquantes qui limite nos conclusions mais aussi remarquer que, malgré cela, de nombreuses perspectives nouvelles de recherche ont émergé qui illustrent la puissance des outils offerts par la génétique des populations.

Conclusions sur les borrélies et I. ricinus en Suisse : 2^e édition

215Une première observation générale à faire ici est que les nouvelles analyses effectuées à l’occasion de cette réédition à l’aide des corrections de l’effet des allèles nuls par FreeNA (Chapuis et Estoup, 2007) (voir aussi le paragraphe correspondant dans la 1^re partie de ce manuel) donne des résultats différents, beaucoup plus cohérents et faciles à interpréter que lors des analyses initiales de la 1^re édition (sans corriger pour l’effet des allèles nuls) avec des estimateurs de subdivision (F_ST) biaisés. Il apparaît donc indispensable de corriger pour les effets des allèles nuls avant de faire des inférences sur la dispersion. Il apparaît également que le mélange de cohortes différentes (ici les échantillons de 1995 et 1996) est loin d’être idéal et qu’il convient également de contrôler pour les effets temporels comme nous l’avons fait lors des analyses réalisées pour la réédition de cet ouvrage. Enfin, il est également apparu que de rendre homozygotes chez les mâles les loci liés au chromosome X (analyses de la 1^re édition) n’était pas une bonne idée et qu’il valait mieux ne prendre en compte les loci hétérosomaux que chez le sexe homogamétique ou supprimer ces loci.

216Quand nous corrigeons pour l’effet des allèles nuls, nous observons un clair isolement par la distance sans avoir recours à des méthodes sophistiquées et complexes de clustérisations bayésiennes. Avec les nouvelles méthodes d’estimation d’effectifs efficaces, nous confirmons l’existence de fortes densités efficaces d’I. ricinus en Suisse (de 600 à 4 000 tiques par km²). Cela conduit à conclure à une très importante viscosité de l’environnement, avec des distances de dispersion de moins de 300 m par génération (donc tous les deux ou trois ans).

217Même en considérant une distribution hétérogène des tiques en Suisse, avec un boisement à 32 %, et une densité par bois moyenne de 859 tiques/km² (celle du bois du Staadswaald), on obtient une densité de 275 tiques/km² sur toute la Suisse (41 285 km²). Les distances de dispersion correspondantes ne dépassent alors pas le km par génération.

218Il existe bien un biais de structuration sexe spécifique. Eu égard à l’ensemble des résultats observés, il est probable que ce biais vienne d’une survie plus faible des femelles immigrantes par rapport aux mâles. Les femelles immatures semblent se nourrir plus souvent chez des hôtes porteurs de B. burdorferi ss, qui en Suisse serait l’écureuil, un hôte peu dispersant. Mais plusieurs travaux ont suggéré la présence de cette borrélie sur d’autres types d’hôtes, en particulier les oiseaux. Il est cependant possible qu’en Suisse l’écureuil soit l’hôte le plus susceptible à cette borrélie et que les femelles immatures d’I. ricinus préfèrent cet hôte. L’écureuil roux disperse rarement à plusieurs km et restreint ses déplacements dans la limite du km (Hämäläinen et al., 2018). Cela pourrait contribuer en partie au biais de structuration observé chez les femelles de la tique. Cependant, le fait que les tiques infectées par cette borrélie montrent très peu de différenciation génétique ne confirme pas un rôle majeur de cette interprétation.

219Il existe globalement une différenciation entre tiques infectées par une ou des borrélies et les tiques saines, ce qui suggère une raciation par l’hôte comme cela a été démontré par Kempf et al. (2011).

220Pour B. afzelii et B. burdorferi non déterminée, il existe une claire signature de biais de structuration pathogène spécifique. Comme les indéterminées sont probablement des borrélies d’oiseaux, une association préférentielle des tiques infectées par ces espèces avec leur hôte vertébré (manipulation de la tique par la bactérie) semble peu convaincante. Il est plus aisé de comprendre ce résultat si on suppose que la survie des tiques immigrantes infectées par ces deux types de borrélies est moins grande que la survie des autres tiques. Cette explication à le mérite d’être cohérente avec la stabilité et la localisation spatiale (endémicité) apparente des foyers de borréliose en Europe (Zeman et Daniel, 1999 ; Randolph, 2001).

221Les borrélies ont une distribution agrégée dans l’espace : certains sites comprenant de nombreuses tiques infectées par l’une des espèces alors que les autres sites en ont très peu. Les différentes espèces ont par ailleurs tendance à s’exclure. Cela provient peut-être d’une interaction entre contraintes spatiales (survie moindre des tiques migrantes et infectées ou transmission moins efficace des tiques immigrantes, présence/absence de certains hôtes) et historiques des différents sites étudiés.

222Si, à la faible efficacité de dispersion de ces borrélies (ou des tiques qui les transportent), d’un site à l’autre, nous ajoutons la forte viscosité de l’environnement pour les tiques que nous avons rappelée plus haut, la relative stabilité et la distribution mosaïque des foyers de borréliose s’expliquent assez facilement.

223La conclusion reste que la biologie des populations de cette espèce de tique et des borrélies qu’elle transmet est très complexe, car dépends de nombreux facteurs en interaction : géographie, sexe de la tique, espèces hôtes rencontrées, infection ou non par tel ou tel autre type de borrélie, etc.

224Une autre conclusion est que nous avons besoin d’affiner ces résultats avec plus de marqueurs de meilleure qualité (avec moins de problèmes d’amplification, autosomaux), ainsi que des échantillons de tiques sur leurs hôtes, ou à tout le moins un barcoding efficace des repas de sang passés, ainsi que, si c’est possible, un outil de détermination du sexe des larves et des nymphes. Cela permettrait en effet de tester nos différentes interprétations avec beaucoup plus d’efficacité.