Informations

Existe-t-il une mesure du degré de réticulation dans un réseau phylogénétique ?


Lorsque l'on travaille avec des réseaux phylogénétiques, comment parler de branches réticulées (c'est-à-dire de type Web) ? Existe-t-il une mesure standard de cela?


Généalogies : les généalogies et les phylogénies sont des réseaux réticulés et pas seulement des arbres divergents

Les pedigrees illustrent les relations généalogiques entre les individus, et les phylogénies font de même pour les groupes d'organismes (comme les espèces, les genres, etc.). Ici, je donne un bref aperçu des concepts et méthodes actuels de calcul et d'affichage des relations généalogiques. Ces relations sont reconnues depuis longtemps comme étant réticulaires, plutôt que strictement divergentes, de sorte que les pedigrees et les phylogénies sont correctement traités comme des réseaux plutôt que comme des arbres. Cependant, actuellement, la plupart des pedigrees sont plutôt présentés comme des « arbres généalogiques » et la plupart des phylogénies sont présentées comme des arbres phylogénétiques. Néanmoins, le développement historique des concepts montre que les réseaux sont antérieurs aux arbres dans la plupart des domaines de la biologie, y compris l'étude des pedigrees, la théorie de la biologie et la pratique de la biologie, ainsi qu'en linguistique historique dans les sciences sociales. Les arbres ont en fait été introduits afin de fournir un modèle conceptuel plus simple pour les relations historiques, puisque les arbres sont un type spécifique de réseau simple. Du point de vue informatique, les arbres et les réseaux font partie de la théorie des graphes, constitués de nœuds connectés par des arêtes. Dans ce contexte mathématique, ils diffèrent uniquement par l'absence ou la présence de nœuds de réticulation, respectivement. Il existe deux types de graphes qui peuvent être appelés réseaux phylogénétiques : (1) les réseaux évolutifs enracinés et (2) les réseaux d'affinité non enracinés. Il existe de nombreuses méthodes de calcul pour les réseaux non enracinés, qui ont deux rôles principaux en phylogénétique : (a) elles agissent comme une forme générique d'affichage de données multivariées et (b) elles sont utilisées spécifiquement pour représenter les réseaux d'haplotypes. Les réseaux évolutionnaires sont plus difficiles à déduire et à analyser, car il n'existe pas d'algorithme mathématique pour reconstruire des événements historiques uniques. Il n'existe donc actuellement aucun cadre d'analyse cohérent pour le calcul de tels réseaux.

Ceci est un aperçu du contenu de l'abonnement, accessible via votre institution.


Préliminaires

Étant donné un ensemble de taxons X . Un sous-ensemble de X (sauf à la fois ∅ et X ) est appelé un cluster. Pour deux grappes C1 et C2 sur X , s'ils sont disjoints ou que l'un est un sous-ensemble de l'autre, c'est-à-dire C1C2 = ∅ ou C1C2 ou C2C1, on dit que C1 et C2 sont compatibles, sinon ils sont incompatibles. Étant donné un ensemble de clusters C sur X . Si deux clusters en C sont compatibles, C est appelé compatible, sinon il est incompatible. Le graphique d'incompatibilité IG( C ) = (V, E) de C est défini comme un graphe non orienté avec un ensemble de nœuds V = C et ensemble de bords E, où deux clusters sont connectés par une arête si et seulement s'ils sont incompatibles.

Laisser S être un sous-ensemble de X . La restriction de C à S, noté C | S , est défini comme le résultat de la suppression de tous les taxons de X S (c'est-à-dire les taxons de X mais pas de S) de chaque cluster dans C . Le degré d' incompatibilité de C , noté d ( C ) , est défini comme le nombre d' arêtes dans I G ( C ) . Soit C un ensemble de clusters sur X et X un taxon dans X . Le degré d'incompatibilité de X, noté par (X), est défini comme le résultat de la soustraction du degré d'incompatibilité de C | X < x > de celui de C , c'est-à-dire d ( x ) = d ( C ) - d ( C | X < x >) . Si le degré d'incompatibilité d'un taxon X est maximale parmi tous les taxons de X , c'est-à-dire d ( x ) = max < d ( y ) | y ∈ X >, nous appelons cela X est le taxon d'incompatibilité w.r.t. C. La fréquence d'un taxon X w.r.t. C , noté F(X), est défini comme le nombre de clusters qui contiennent X, c'est-à-dire f ( x ) = | < C ∈ C | x ∈ C >| .

Un réseau phylogénétique enraciné N = (V, E) sur X est un graphe acyclique orienté enraciné (DAG en abrégé), et ses feuilles sont étiquetées bijectivement comme X . Le degré d'un nœud vV est noté comme σ(v). Un nœud v est un nœud réticulé si σ(v) ≥ 2 sinon c'est un nœud d'arbre en particulier, un nœud d'arbre est un nœud racine si σ(v) = 0. Une arête e = (vous, v) est une arête d'arbre si v est un nœud d'arbre, sinon c'est une arête réticulée. Le numéro de réticulation dans un réseau N = (V, E) est ∑σ(v(σ(v) − 1) = |E| − |V| + 1.

Étant donné un cluster C et un arbre phylogénétique enraciné T. S'il y a un bord e dans T tel que l'ensemble des taxons accessible à partir de e équivaut à C, on dit que T représente C. Étant donné un réseau N, lors de la connexion d'un bord entrant et de la déconnexion de tous les autres bords entrants pour chaque nœud réticulé, s'il existe un bord d'arbre e tel que l'ensemble de taxons accessible à partir de e équivaut à C, on dit que N représente C au sens câblé. Sinon, s'il y a une lisière d'arbre e dans N tel que l'ensemble de taxons accessible à partir de e équivaut à C, on dit que N représente C au sens câblé.

Laisser N = (V, E) être un réseau représentant l'ensemble des clusters C . Un cluster C ∈ C est souvent représenté par plus d'une arête d'arbre dans N et une lisière d'arbre eE représente souvent plus d'un cluster dans C . S'il existe un mappage ϵ de C à l'ensemble des arêtes de l'arbre de N, où ϵ(C) est une arête d'arbre représentant C pour C ∈ C , tel que pour deux clusters C 1 , C 2 ∈ C , C1 et C2 se trouvent dans la même composante connexe du graphe d'incompatibilité IG( C ) si et seulement si ϵ(C1) et ϵ(C2) sont contenus dans le même composant biconnecté de N. Alors on appelle ça N est décomposable w.r.t. C.

Lors de la construction de réseaux phylogénétiques enracinés à partir d'arbres phylogénétiques enracinés, les procédés calculent d'abord les clusters représentés par les arbres d'entrée, puis construisent un réseau phylogénétique enraciné représentant tous les clusters.

Les réseaux phylogénétiques enracinés peuvent décrire l'histoire de l'évolution en présence d'événements réticulés, tels que les transferts horizontaux de gènes, les hybridations et les recombinaisons. Ces événements réticulés sont rares en réalité [33]. En conséquence, on s'attend à ce que le réseau construit ait le nombre minimal de nœuds réticulés. Laisser N être un réseau construit pour l'ensemble de clusters d'entrée C . Supposons que C ′ est l'ensemble des clusters représentés par N. En fait C ′ contient plus de clusters que C , c'est-à-dire C ⊊ C ′ . On définit ici les clusters redondants C 0 de N comme les clusters qui sont dans C ′ mais pas dans C , c'est-à-dire C 0 = C ′ C . En analyse phylogénétique, les taxons d'un groupe sont présumés monophylétiques. Par conséquent, la situation idéale serait C 0 = ∅ , c'est-à-dire que tous les clusters représentés par le réseau construit seraient les clusters représentés dans les arbres d'entrée, et aucun autre. Par conséquent, au moyen du principe de parcimonie, le réseau le mieux construit est celui qui minimise le nombre de clusters redondants, ce qui est basé sur la condition préalable qu'il ait un nombre minimal de nœuds réticulés.

La section suivante présentera les principales méthodes de construction de réseaux phylogénétiques enracinés à partir d'arbres phylogénétiques enracinés.


Résultats

Densité de probabilité d'un arbre génétique.

Étant donné un réseau phylogénétique Ψ et une généalogie génétique G j pour le locus j (topologie et longueurs de branches dans les deux cas), on note H Ψ ( G j ) l'ensemble de toutes les histoires coalescentes de G j au sein des branches de Ψ. Ensuite, la distribution (densité) des arbres de gènes est donnée par p ( G j | Ψ , Γ ) = ∑ h ∈ H Ψ ( G j ) p ( h | Ψ , Γ ) , [4] où Γ est les probabilités d'héritage matrice, comme décrit ci-dessus. Pour une arête b = ( x , y ) ∈ E ( Ψ ) , nous définissons T b ( h ) comme étant le vecteur des temps (par ordre croissant) des événements de coalescence qui se produisent sur la branche b sous l'histoire coalescente h et l'heure du nœud oui (une définition formelle est fournie dans Annexe SI). On note T b ( h ) [ i ] le jeème élément du vecteur. De plus, on note u b ( h ) le nombre de lignées génétiques entrant dans le bord b et par v b ( h ) le nombre de lignées génétiques quittant le bord b sous h. Alors, on a p ( h | Ψ , Γ ) = ∏ b ∈ E ( Ψ ) [ ∏ i = 1 | T b ( h ) | − 1 e − ( ub ( h ) − i + 1 2 ) ( T b ( h ) i + 1 − T b ( h ) i ) ] × e − ( vb ( h ) 2 ) ( τ Ψ ( b ) − T b ( h ) | T b ( h ) | ) × Γ [ b , j ] ub ( h ) , [5] où τ Ψ ( b ) pour l'arête b = ( x , y ) est le temps du nœud X dans le réseau phylogénétique Ψ. Une dérivation complète de la formule et un algorithme plus efficace pour le calculer le long des lignes de Yu et al. (17), qui évitent des sommations explicites sur les histoires coalescentes possibles, sont donnés dans Annexe SI.

Recherche dans l'espace des réseaux phylogénétiques.

Soit Ω ( n ) l'espace de tous les réseaux phylogénétiques sur m taxons, on note Ω ( n , k ) le sous-espace de Ω ( n ) qui contient tous les réseaux phylogénétiques (rDAG) avec m feuilles et k nœuds de réticulation. En particulier, Ω ( n , 0 ) est le sous-espace qui contient tous les arbres phylogénétiques. Pour rechercher l'espace du réseau phylogénétique en couches, nous définissons deux opérations qui permettent de rechercher dans Ω ( n , k ) un k: une opération qui permet à la recherche de remonter une couche de ( n , k ) à Ω ( n , k + 1 ) et une opération qui permet à la recherche de descendre une couche de ( n , k ) à Ω ( n , k−1 ). Pour la recherche au sein d'une couche, les opérations déplacent soit la destination d'une arête de réticulation, soit la source d'une arête (réticulation ou non). Pour remonter une couche, l'opération consiste à ajouter une arête de réticulation entre deux arêtes existantes dans le réseau, et pour descendre une couche, l'opération supprime une arête de réticulation (plus de détails sont fournis dans Annexe SI). Il convient de mentionner que bien que l'espace de toutes les topologies d'arbres phylogénétiques sur m taxons est fini, l'espace de toutes les topologies de réseaux phylogénétiques sur m taxons est, en théorie, infini, car Ω ( n ) = ∪ k ≥ 0 Ω ( n , k ) et k sont illimités. Par exemple, considérons le cas de seulement deux taxons. Il y a un arbre enraciné unique dans ce cas. Cependant, étant donné que plusieurs événements d'hybridation peuvent se produire entre les deux mêmes taxons frères à des moments différents, n'importe quel nombre de bords horizontaux peut être ajouté entre ces deux taxons. Néanmoins, la question de savoir si de tels scénarios d'hybridation répétitive sont identifiables à partir d'ensembles de données génomiques typiques est une autre question.

Une heuristique pour estimer les longueurs de branches optimales pour une topologie d'arbre d'espèces fixe, étant donné les topologies d'arbres de gènes, qui est basée sur l'application répétée de la méthode de Brent (22) a été introduite par Wu (20). Nous utilisons une heuristique similaire pour estimer les longueurs de branches du réseau phylogénétique et les probabilités d'héritage (tous les détails sont donnés dans Annexe SI). Le couplage des transformations topologiques et des heuristiques d'estimation des paramètres avec la formulation de vraisemblance ci-dessus permet de rechercher l'espace de manière ascendante pour déduire un réseau phylogénétique ML. Étant donné l'existence d'optimums locaux au sein de chaque couche, de multiples analyses indépendantes peuvent être effectuées.

Contrôle de la complexité du modèle.

Parce que les réseaux dans ( n , k + 1 ) fournissent des modèles plus complexes que les réseaux dans ( n , k ) , les approches décrites ci-dessus doivent gérer le problème de sélection de modèle. Des critères d'information ont déjà été utilisés dans le cadre de réseaux phylogénétiques (12, 16), et nous les utilisons ici (au lieu de rechercher sur la base du score de vraisemblance, la recherche procède sur la base des valeurs de ces critères, qui intègrent les scores de vraisemblance) . Une autre approche que nous proposons ici, pour la première fois à notre connaissance, est l'utilisation de K-fold validation croisée, par laquelle l'ensemble d'arbres génétiques d'entrée est partitionné en K sous-ensembles de tailles égales, les paramètres du modèle sont déduits de K − 1 sous-ensembles, et l'ajustement du modèle du sous-ensemble restant est calculé. Cet ajustement est calculé en comparant les fréquences des arbres de gènes dans le sous-ensemble de validation avec la distribution des arbres de gènes produits par le réseau inféré. Si l'ajustement du meilleur réseau Ψ ″ trouvé dans Ω ( n , k + 1 ) n'est pas bien meilleur (nous utilisons un seuil d'amélioration de 3%, choisi sur la base d'observations empiriques) que l'ajustement du meilleur réseau Ψ ′ trouvé dans Ω ( n , k ) , on déclare k être l'estimation correcte du nombre de nœuds de réticulation et Ψ ′ être le réseau phylogénétique optimal. Il est important de noter ici que cette idée de validation croisée ne fonctionne que pour les topologies d'arbres de gènes entièrement résolues, car dans le cas d'arbres de gènes avec des longueurs de branches, les fréquences des arbres de gènes dans le sous-ensemble de validation ne sont pas informatives.

Enfin, pour évaluer le support des réseaux phylogénétiques que nous inférons, nous proposons d'utiliser le bootstrap paramétrique. Après avoir inféré un réseau Ψ à partir des données G , nous utilisons Ψ pour générer ℓ jeux de données, à partir desquels nous déduisons ℓ réseaux phylogénétiques Ψ 1 , … , Ψ ℓ . Nous estimons ensuite le support de chaque branche b dans Ψ comme le nombre de réseaux (sur les ℓ ) qui ont une branche équivalente. On dit que deux arêtes dans deux réseaux phylogénétiques sont équivalentes si (je) l'un ou les deux sont des bords de réticulation ou les deux ne sont pas et (ii) définissent tous deux les mêmes groupes (le groupe défini par une branche est l'ensemble de tous les taxons de cette branche dans le réseau).

Performances sur données simulées.

Nous avons mis en œuvre toutes les méthodes décrites ci-dessus dans le progiciel open source disponible au public PhyloNet (23) et étudié les performances des méthodes sur plusieurs ensembles de données simulés. Dans l'étude de simulation dont les résultats sont rapportés sur la figure 2, nous avons utilisé le réseau phylogénétique 1 comme réseau modèle, et pour différents nombres de loci, nous avons fait évoluer des arbres de gènes sous la coalescence au sein des branches du réseau, puis simulé l'évolution des séquences sur ces arbres de gènes avec des longueurs de séquences variables. Nous avons ensuite estimé pour chaque alignement de séquences 100 arbres de gènes en utilisant ML avec bootstrapping. Enfin, nous avons déduit des réseaux en utilisant notre méthode ML à partir de (je) de véritables topologies d'arbres de gènes, (ii) les topologies estimées des arbres de gènes, (iii) les vraies topologies d'arbres génétiques et les longueurs de branches, et (iv) a estimé les topologies d'arbres génétiques et les longueurs de branches. Les résultats de (je) et (ii) sont illustrés à la figure 2B, alors que les résultats de (iii) et (iv) sont illustrés à la figure 2C. Pour chaque réglage du nombre de loci et de la longueur de séquence, nous avons généré 30 ensembles de données et effectué des inférences sur chacun d'eux.

Précision de la méthode sur les données simulées. (UNE) Les données ont été générées le long du réseau phylogénétique Ψ 1 (toutes les branches internes, à l'exception du bord horizontal, ont des longueurs de 1 unité coalescente et la probabilité d'héritage est de 0,1 pour tous les loci). Résultats basés sur des estimations de la topologie de l'arbre génétique (B) et les estimations de la topologie de l'arbre génétique et de la longueur des branches (C) sont indiqués. Pour chaque nombre de loci, la barre la plus à droite correspond à l'inférence à partir des véritables généalogies génétiques et les trois autres barres, de gauche à droite, correspondent aux généalogies génétiques estimées (en utilisant 100 répliques bootstrap et Eq. 3) à partir de séquences de longueurs de 250, 500 et 1 000, respectivement. Les régions bleu foncé, cyan et jaune correspondent au nombre de fois que chacun des réseaux Ψ 1 , 2 et Ψ 3 , respectivement, dans UNE a été déduit. La région marron correspond au nombre de fois qu'un autre réseau avec une seule réticulation a été déduit. Ici, un individu a été échantillonné par taxon pour chacun des loci.

Alors que l'hybridation dans le réseau modèle implique B et l'ancêtre commun le plus récent (MRCA) de C et D, la longueur de la branche entre l'événement d'hybridation et la divergence de C et D de leur MRCA peut avoir un effet important sur la distinction entre les vrais scénarios d'hybridation et les deux donnés par Ψ 2 et Ψ 3 dans la Fig. 2UNE. Par conséquent, pour chaque ensemble de données, nous avons enregistré si la méthode inférait l'un des trois réseaux illustrés à la Fig. 2UNE, contrairement à tout autre réseau à réticulation unique.

Plusieurs tendances peuvent être observées sur la figure 2UNE. Premièrement, l'utilisation des véritables topologies d'arbres de gènes avec des longueurs de branches permet d'obtenir des inférences plus précises que l'utilisation de topologies d'arbres de gènes seules. Ce résultat n'est pas surprenant, car le premier type de données contient plus d'informations que le second. En particulier, lors de l'utilisation de 80 ou 160 loci, le réseau déduit des vrais arbres de gènes avec des longueurs de branches est toujours le vrai réseau. D'autre part, en utilisant uniquement les topologies d'arbres génétiques pour 160 loci, dans cinq des 30 cas, l'inférence a renvoyé l'un des deux réseaux alternatifs Ψ 1 et Ψ 2 . Deuxièmement, la précision des inférences s'améliore à mesure que le nombre de loci augmente et que la longueur de la séquence augmente, bien que l'augmentation du nombre de loci ait eu un effet beaucoup plus positif sur la précision de l'inférence. Troisièmement, un résultat très surprenant est qu'en utilisant uniquement des topologies d'arbres de gènes, l'utilisation des vrais arbres de gènes n'a presque jamais donné une meilleure précision que lors de l'utilisation de topologies d'arbres de gènes estimées pour un nombre donné de loci. Ce résultat atteste du fait qu'en tenant soigneusement compte de l'incertitude dans les estimations de l'arbre génétique, la méthode peut obtenir de très bons résultats. Même en utilisant des topologies d'arbres de gènes et des longueurs de branches, le gain de précision lors de l'utilisation des vrais arbres de gènes est très faible par rapport à l'utilisation des estimations d'arbres de gènes avec l'incertitude prise en compte. Quatrièmement, la combinaison d'une faible valeur de probabilité d'hérédité (0,1 dans cette simulation) et d'un temps relativement court entre l'hybridation et la spéciation ultérieure entraîne une incertitude dans l'identification du donneur et du receveur de l'événement d'hybridation. Par exemple, lors de l'utilisation de topologies d'arbres génétiques uniquement pour 160 loci, le réseau inféré est toujours l'un des trois réseaux Ψ 1 , Ψ 2 et Ψ 3 , même s'il s'agit principalement de Ψ 1 . Nous avons constaté que l'augmentation de la longueur des branches ou des probabilités d'héritage se traduirait par des précisions plus élevées. De plus, dans nos simulations, nous avons constaté que l'augmentation du nombre d'individus échantillonnés par taxon entraînerait une amélioration de la précision, quoique plutôt légèrement (Annexe SI). Cependant, nous nous attendons à ce que l'échantillonnage d'un plus grand nombre d'individus se traduise par des améliorations plus significatives sur des ensembles de données plus grands ou plus complexes. En termes de probabilités d'hérédité inférées, les vrais arbres de gènes ont donné des estimations très précises, tandis que les arbres de gènes estimés avec des longueurs de branches ont abouti, en général, à des estimations plus précises des probabilités. Enfin, nous avons constaté que la validation croisée fait généralement mieux que les critères d'information pour déterminer le nombre d'événements d'hybridation (y compris sur l'ensemble de données biologiques, comme discuté ci-dessous). Des résultats de simulation plus complets dans des scénarios plus faciles à inférer que ceux dont nous avons discuté ici sont contenus dans Annexe SI.

Analyse d'un ensemble de données de souris domestiques multilocus.

Nous avons également utilisé notre méthode pour analyser un ensemble de données multilocus de souris domestique (M. musculus) des génomes, obtenus à partir des études de Staubach et al. (7), Didion et al. (24), et Yang et al. (25) (plus de détails sont fournis dans Annexe SI). Staubach et al. (7) ont trouvé des preuves substantielles à l'échelle du génome d'introgression sous-spécifique dans les quatre populations, représentant 3 % du génome dans les deux M. m. domestique (l'une de France et l'autre d'Allemagne), 4 % dans un M. m. musculature population du Kazakhstan, et 18% dans un M. m. musculature population de la République tchèque. Cependant, il est important de noter que la méthode HAPMIX (26), qui a été utilisée par Staubach et al. (7), ne tient pas explicitement compte de l'ILS.

Notre étude a inclus tous les échantillons de l'étude de Staubach et al. (7). De plus, notre étude a inclus des échantillons supplémentaires provenant d'un M. m. musculature population de Chine (25) qui n'ont pas été utilisées dans l'étude de Staubach et al. (7). Dans cette analyse, nous avons utilisé uniquement des topologies estimées d'arbres de gènes. La raison en est que les séquences génomiques sont obtenues à partir d'individus très proches (ces individus sont cinq individus de la même espèce), et très peu de variation existe dans les données pour estimer la longueur des branches avec précision. De plus, cette faible variation ne permet pas une analyse bootstrap appropriée des arbres de gènes pour les loci individuels. Le signal puissant de cet ensemble de données provient du très grand nombre de loci. Dans notre analyse, nous avons trouvé une amélioration significative dans un réseau phylogénétique avec une seule réticulation sur aucune réticulation et une amélioration significative dans un réseau phylogénétique avec deux réticulations sur une seule réticulation. Cependant, lorsque nous avons poursuivi la recherche du réseau optimal avec trois réticulations, nous avons constaté que l'amélioration obtenue en considérant un troisième événement de réticulation était insignifiante sur la base des critères d'information, et qu'il n'y avait aucune amélioration du tout sur la base de la validation croisée. Nous avons ainsi appelé le réseau phylogénétique optimal avec deux réticulations comme hypothèse pour l'histoire évolutive de cet ensemble de génomes. Cette histoire évolutive est illustrée à la Fig. 3 (plus de détails sur les résultats et les analyses sont fournis dans Annexe SI). Le réseau phylogénétique n'est pas ultramétrique, et il convient de souligner que les longueurs de branches sont données en unités coalescentes. Ainsi, le manque d'ultramétrie pourrait être dû à des tailles de population différentes ou, dans une moindre mesure, à des temps de génération différents.

Réseau phylogénétique optimal déduit sur la souris domestique (M. musculus) base de données. Un seul individu a été échantillonné dans chacune des cinq populations : M. m. domestique de France (DF), M. m. domestique d'Allemagne (DG), M. m. musculature de la République tchèque (MZ), M. m. musculature du Kazakhstan (MK), et M. m. musculature de Chine (MC). L'analyse a trouvé plusieurs réseaux phylogénétiques presque également optimaux avec deux événements de réticulation. Ces réseaux multiples étaient tous d'accord sur les populations bénéficiaires mais en désaccord sur les populations donneuses. Une hybridation (la flèche horizontale en pointillés en haut) implique le MRCA de DF et DG en tant que population receveuse, mais semble avoir impliqué MK, MC ou leur MRCA en tant que population donneuse. La seconde hybridation (la flèche horizontale en pointillés du bas) implique MZ en tant que population receveuse, mais semble avoir impliqué DF, DG ou leur MRCA en tant que population donneuse. Les longueurs de branches en unités coalescentes (sur les branches de l'arbre) et les probabilités d'héritage (sur les bords horizontaux) sont indiquées (tous les détails des données et des résultats sont fournis dans Annexe SI).

Notre analyse des génomes des souris domestiques produit une histoire évolutive qui diffère de celle rapportée par Staubach et al. (7) non seulement en termes de nombre de populations impliquées mais aussi en prenant en compte l'histoire évolutive des populations impliquées. Nous considérons les pourcentages du génome d'origine introgressée rapportés par Staubach et al. (7) être des surestimations, car l'introgression impliquant une population ancestrale qui s'est ensuite divisée en plus d'une population existante serait rapportée à plusieurs reprises pour chaque population existante dans le cas de l'étude de Staubach et al. (7). En revanche, les mêmes pourcentages seraient sous-estimés dans le cas où des populations mélangées seraient utilisées à la place des populations de référence non mélangées requises par HAPMIX, comme Staubach et al. (7) a fait en utilisant des échantillons de souris putativement introgressés pour construire les populations de référence. Notamment, notre méthodologie ne nécessite pas l'utilisation de populations de référence non mélangées.

Nous émettons l'hypothèse que l'événement d'introgression le plus récent de la figure 3 est dû au flux de gènes provenant du contact secondaire, où les plages des deux M. musculus sous-espèces se chevauchaient, à peu près à la frontière entre l'Allemagne et la République tchèque. L'interprétation biologique de l'événement d'introgression plus ancien est moins claire. Nous conjecturons que l'événement est lié au flux de gènes pendant et après la divergence sous-spécifique. Une étude plus approfondie peut fournir des indices importants sur la base mécaniste de l'évolution des sous-espèces dans M. musculus et le processus de spéciation lui-même.

Il est important de noter que bien que nous ayons utilisé un très grand nombre de loci, il restait une incertitude quant aux origines inférées des deux événements d'hybridation (comme le montre la figure 3), un schéma similaire à celui observé dans les résultats de la simulation et discuté ci-dessus. Cette incertitude est le reflet du signal faible dans ces données, couplé à la faible probabilité d'hérédité et à la courte longueur de branche entre l'hybridation et le MRCA de M. m. musculature de Chine et M. m. musculature du Kazakhstan et M. m. domestique de France et M. m. domestique d'Allemagne, ce qui est une question que nous avons abordée ci-dessus dans le contexte des données simulées. Les échantillons utilisés sont très proches, ce qui donne des génomes avec un très petit nombre de sites de ségrégation, et donc un signal plus faible pour l'inférence. Néanmoins, l'incertitude est localisée dans le sens où les donneurs potentiels du matériel génétique de chaque événement d'hybridation tournent autour d'un seul nœud ancestral. Étant donné que les cinq populations analysées sont étroitement liées, la plupart des arbres génétiques reconstruits n'étaient pas binaires, en raison de séquences identiques d'allèles multiples. Parce que l'amorçage n'est pas utile dans ces scénarios (chaque locus a une poignée de sites, dont la plupart sont monomorphes), nous avons utilisé les topologies d'arbres de gènes non binaires pour les loci et considéré l'ensemble de toutes les résolutions comme l'ensemble d'estimations d'arbres de gènes à utiliser dans l'éq. 3.


Fond

Les réseaux phylogénétiques sont des graphiques utilisés pour représenter les relations phylogénétiques entre différents taxons, et sont généralement utilisés lorsqu'une représentation arborescente ne suffit pas. Il existe de nombreux types de réseaux phylogénétiques et il est utile de distinguer deux classes principales : implicite les réseaux phylogénétiques qui fournissent des outils pour visualiser et analyser les signaux phylogénétiques incompatibles, tels que les réseaux divisés [1, 2] et explicite les réseaux phylogénétiques qui fournissent des scénarios explicites d'évolution réticulée, tels que les réseaux d'hybridation [3–7], les réseaux HGT [8] et les réseaux de recombinaison [9–19].

Le logiciel actuellement disponible pour le calcul et l'analyse de réseaux phylogénétiques explicites consiste en une diffusion d'implémentations de base d'algorithmes développés pour résoudre la tâche de calcul [6, 14, 16-18, 20]. La plupart des logiciels sont pilotés par ligne de commande et une visualisation attrayante des résultats fait souvent défaut. Il est essentiel de disposer d'un outil qui permet à la fois une large utilisation des méthodes disponibles pour les biologistes, et un développement meilleur et plus poussé de nouvelles méthodes.

SplitsTree est une application développée dans notre groupe de recherche, visant à l'origine l'analyse phylogénétique d'ensembles de scissions. La dernière version de SplitsTree [2] intègre une variété de méthodes pour le calcul, la visualisation et l'interprétation des arbres phylogénétiques et des réseaux phylogénétiques implicites. Deux principaux avantages de SplitsTree sont l'interface utilisateur graphique (GUI) et l'intégration d'algorithmes via un chargeur de classe piloté par interface (plugins). Dans cet article, nous présentons une extension de SplitsTree qui permet au programme de gérer des réseaux phylogénétiques explicites. L'extension résout deux problèmes importants : une intégration efficace des réseaux phylogénétiques explicites, et la visualisation de ces réseaux.


Méthodes

Liu et al. a récemment introduit MP-EST, une approche de pseudo-vraisemblance maximale pour estimer les arbres d'espèces à partir d'une collection d'arbres de gènes enracinés sous la fusion multi-espèces [26]. La méthode a permis d'améliorer considérablement le temps d'exécution de l'inférence statistique des arbres d'espèces. Inspiré par ce travail, nous proposons une méthode pour estimer la phylogénie des espèces en présence à la fois d'hybridation et de tri de lignée incomplet sous une pseudo-vraisemblance maximale.

Réseaux phylogénétiques, arbres de gènes et triplets enracinés

Un réseau phylogénétique (binaire) [29] Ψ sur l'ensemble X d'espèces (taxons) est un graphe enraciné, orienté, acyclique dont l'ensemble de nœuds est V (Ψ) = <r> V L V T V N

r est la racine de et satisfait (r) = 0 et + (r) = 2

V L: l'ensemble de feuilles étiqueté bijectivement par X , où (v) = 1 et + (v) = 0 pour tout v V L

V T: nœuds internes de l'arbre, où (v) = 1 et + (v) = 2 pour tout v VT et,

V N: nœuds de réticulation, où (v) = 2 et + (v) = 1 pour tout v V N.

Ici, (v) et + (v) sont les degrés d'entrée et de sortie du nœud v, respectivement. On note par E(Ψ) l'ensemble des arêtes du réseau Ψ. Le réseau phylogénétique a des longueurs de branches ?? : E(Ψ) + . Ci-après, nous utiliserons pour désigner à la fois la topologie et la longueur des branches d'un réseau phylogénétique. De plus, comme dans [22, 24], pour un cadre probabiliste, il existe une fonction supplémentaire, appelée probabilité d'héritage, ?? : E(Ψ) [0, 1] qui satisfait :

??(e) = 1 pour chaque arête e dont la tête est un nœud d'arbre, et

??(e1) + ??(e2) = 1 pour chaque paire d'arêtes e1 et e2 dont la tête est le même nœud de réticulation.

Dans [24], nous avons discuté de la façon de généraliser la fonction ?? de sorte qu'il varie selon les loci, et cette généralisation serait triviale à incorporer dans les méthodes ci-dessous.

Un arbre génétique g sur l'ensemble X d'espèces est un arbre enraciné (pas nécessairement binaire) dont les feuilles sont étiquetées (pas nécessairement de manière bijective) par X . Pour distinguer les feuilles étiquetées par le même élément de X , nous ajoutons des indices aux étiquettes des feuilles. La figure 1 montre un arbre génétique g sur le plateau X = <X, Oui, Z> d'espèces, où quatre allèles sont échantillonnés à partir d'espèces X (étiqueté X1, . X4), trois allèles sont échantillonnés à partir d'espèces Oui (étiqueté oui1, . oui3), et deux allèles sont échantillonnés à partir d'espèces Z (étiqueté z1 et z2). En particulier, dans ce travail, nous permettons à un arbre génétique d'avoir zéro allèle échantillonné de certaines espèces.

Arbres génétiques et triplets enracinés. Un arbre génétique g sur trois espèces X, Y , et Z, où plusieurs allèles sont échantillonnés par espèce. Deux triplets induits par l'arbre génétique et leur correspondance avec les noms d'espèces sont montrés.

Un triple enraciné (à partir de maintenant on écrira simplement « triple », puisque nous ne traitons que des topologies enracinées) est un arbre enraciné à trois feuilles. Si le triplet est binaire, on écrit xy|z pour indiquer que le triple met X et oui plus proches l'un de l'autre que l'un d'eux pour z. Si le triplet est non binaire, alors il est xyz. On note par g|<x,y,z> le triple dans l'arbre des gènes g induite en restreignant sa foliation aux trois feuilles étiquetées x, y, et z. La figure 1 montre les deux triplets induits par <X1, oui1, z1> et <X1, oui1, z2>. Enfin, pour lier les étiquettes des feuilles dans l'arbre génétique à leurs taxons correspondants dans le réseau phylogénétique, nous introduisons la fonction ?? qui mappe une étiquette d'allèle dans l'arbre génétique à son taxon correspondant dans le réseau. Par exemple, dans la figure 1, ??(z1) = ??(z2) = Z. De plus, nous utilisons ??(g|<x,y,z>) pour désigner le triple induit avec ses étiquettes foliaires remplacées par les noms de taxons (d'espèces). La figure 1 illustre ??.

Pseudo-vraisemblance d'un réseau d'espèces

Soit X un ensemble de taxons (espèces), t = XY |Z être un binaire tripler avec X, Y, Z X , et g être un arbre de gènes sur X . On note par une(g, X), pour XX, l'ensemble des allèles de X cette étiquette laisse de g. Par exemple, dans la figure 1, une(g, X) = <X1, X2, X3, X4>. Nous définissons ??(t, g) pour être le nombre de fois t est induit par g (lorsque les étiquettes de feuilles sont mappées sur X à l'aide de la fonction ??) normalisé par le nombre de fois qu'un triplet sur X, Y , et Z est induit par g. Il est clair que si au plus un allèle par espèce est échantillonné dans g, alors tout triple n'est soit pas induit par l'arbre génétique, soit induit une fois. Cependant, étant donné que nous autorisons plusieurs allèles par espèce, cela pourrait ne pas être le cas. Notez que tandis que t est binaire, il se peut que g|<X je , y j ,z k> est non binaire. Puisqu'il existe trois façons de résoudre un triplet non binaire, un triplet non binaire g|<X je , y j ,z k> contribue 1/3 à la valeur de ??(t, g). Compte tenu de ces deux problèmes, ??(t, g) pour t = XY |Z équivaut à

où I est la fonction indicatrice définie par I e = 1 lorsque e est vrai et I e = 0 quand e c'est faux. Pour un ensemble g des arbres de gènes, nous définissons ρ ( t , G ) = ∑ g ∈ G ρ ( t , g ) . Si le dénominateur dans l'Eq. (1) est égal à zéro, on pose ??(t, g) = 0.

Étant donné un ensemble g des arbres génétiques, les trois triplets binaires t1 = XY|Z, t2 = XZ|Y, et t3 = YZ|X sur un ensemble < X , Y , Z >⊆ X ont une distribution multinomiale donnée par

P (t|??,) est la probabilité d'un triplet enraciné t réseau donné Ψ et probabilités d'héritage ?? [22, 21].

Finally, the pseudo-likelihood of phylogenetic network Ψ and inheritance probabilities ?? given a set g of gene trees is given by

A maximum pseudo-likelihood approach seeks Ψ et ?? that maximize Eq. (3).

Since for a given set g of gene trees | G | ! ∏ i = 1 3 ρ ( t i , G ) ! is a constant irrespective of Ψ and ??, this term is dropped from the pseudo-likelihood computation when searching for Ψ et ?? .

Convergence and identifiability

It follows from the strong law of large numbers [30] that as the number of gene trees

|G| goes to infinity, the proportions of rooted triples in gene trees converge to their

where Ψ ^ is the true phylogenetic network and γ ^ are the true inheritance probabilities. Thus, as |G| goes to infinity, L, γ|G) converges to

A phylogenetic tree is uniquely encoded by its triple system [31]. More specifically, given a phylogenetic tree T , laisser R(T) be the set of triples induced by tree T . Then no tree T′ exists such that T ≠ T′ et R(T) = R(T′). Combining this fact with Eq. (5) and the fact that H,) is maximized when λ = λ ^ and γ = γ ^ , it is clear that when the species phylogeny Ψ is a tree, as |G| goes to infinity, Ψ converges to the true species tree [26].

However, in contrast to trees, triples do not necessarily uniquely encode a phylogenetic network [32]. For example, the three phylogenetic networks Ψ1, Ψ2 et3 in Figure 2 have different topologies, but they induce (a network induces a triple if at least one of the trees displayed by the network induces that triple) the same triple system <A|BC, AB|C, A|BD, AB|D, A|CD, B|CD>. This means that, given a phylogenetic network Ψ (topology and branch lengths) and inheritance probabilities ??, if there is a phylogenetic network Ψ′ s.t. R(Ψ) = R(Ψ′) (which is not necessarily always true), then there exist branch lengths for Ψ′ and inheritance probabilities ??′ such that P (t|??,) = P (t|Ψ′,′) for every rooted triple t. For example, in Figure 2, given network Ψ1 with its branch lengths and inheritance probabilities, we can obtain P (t|??1,1) = P (t|??2,2) for every triple t by setting the branch lengths of network Ψ2 and inheritance probabilities as

Illustration of the lack of network identifiability under the proposed pseudo-likelihood framework. Three phylogenetic networks with the same set of triples: A|BC, AB|C, A|BD, AB|D, A|CD, et B|CD. Branch lengths and inheritance probabilities are shown in blue and red, respectively, for Ψ1 et2.

A concrete example of these settings is:

network Ψ1: b1 = 1, b2 = 1, b3 = 2, b4 = 1, b5 = 0, ?? = 0.1

network Ψ2: je1 = 1.841435, je2 = 1.951019, je3 = 0.207841, ?? = 0.6631633.

This result means that when a species network Ψ is not uniquely encoded by its triple system, as the number of gene trees |G| goes to infinity, argmaxΨ,γ L, γ|G) is not unique, and one of the solutions is the true species network Ψ ^ and true inheritance probabilities γ ^ . This leads to an issue in our inference: if the optimal phylogenetic network Ψ ^ is not uniquely encoded by its triple system R ( Ψ ^ ) , the maximum pseudo-likelihood search might return any of the optimal networks with the same triple system. To ameliorate (yet, not guaranteed to always solve) the identifiability issue, one heuristic is to save all optimal networks identified during the search based on pseudo-likelihood and then optimize their branch lengths and inheritance probabilities using the full likelihood computation [22, 21] to identify the optimal one among them. However, it is important to keep in mind that full likelihood computation can be infeasible except for very small data sets.

Searching for Ψ ∗ et ?? ∗

Given a set of gene trees g, Ψ et ?? that maximize L, γ|G) are searched by traversing the space of phylogenetic networks and inheritance probabilities using simulated annealing. The search starts from initial values of Ψ and ?? and in every iteration, one of the following operations is selected randomly according to their preset weights:

Modifying one or more branch lengths.

Modifying one or more inheritance probabilities.

Adding a reticulation edge.

Deleting a reticulation edge.

Relocating the head of a reticulation edge.

Relocating the tail of an edge.

The first two operations do not change the topology of the network. Full details of how these operations are implemented are given in [24]. During the search, if the new network has higher pseudo-likelihood than the current one, it is always accepted otherwise, it is accepted with some probability. The search terminates if one of two conditions is satisfied: (1) the number of iterations reaches some preset maximum value or (2) the search is alternating between a collection of species networks with high pseudo-likelihoods, and a sufficient number of iterations have passed since visiting any other species networks. The details of how the probability of acceptance is set and the termination conditions are determined are similar to those used in [33, 34].

Since branch lengths and inheritance probabilities are sampled, rather than optimized, during the search, some solutions could be missed due to this sampling. One heuristic to ameliorate this problem is to keep the top k optimal networks during the search, and then at the end optimize the branch lengths and inheritance probabilities (under the pseudo-likelihood criterion) of only these networks to identify the optimal one. We implemented this in our method and we discuss its performance in the simulation results below.


Key Processes of Divergence and Reticulation in Nature

Comparative Phylogeography Across the Australian Monsoonal Tropics.

In essence, comparative phylogeography is about establishing commonalities of spatial patterns of genetic and gene tree diversity across codistributed species (47, 48). Combined with population genetic (coalescent) and spatial modeling (49), this effort has yielded insights into biogeographic history, such as locations of refugia and expansion areas (50) and the varying effects of ecological or physical dispersal barriers. In a comparative setting, such studies can identify how landscape features and regional climatic variation have interacted with the varying ecologies of species to shape current diversity (51) and how these interactions can influence speciation processes (52, 43).

To explore divergence vs. reticulation processes in a comparative context, we focus on phylogeographic data for ecologically diverse species from northern Australia, a vast stretch of monsoonal savanna and woodlands with interspersed ancient sandstone plateaus (54) (Fig. 3). That this rich tropical fauna is biogeographically structured has long been known (55) and formalized using cladistic biogeography by Cracraft (56), who recognized a basal dichotomy across the treeless “Carpentarian Barrier” (CB) separating the Kimberley (KIM) and Top End (TE) faunas from the faunas in Cape York (CY) and the eastern Australian forest (EF), as well as New Guinea. KIM and TE are, in turn, separated by hot, low-relief, and relatively dry regions associated with several smaller barriers (57), collectively referred to here as the Kimberley-Top End Barrier (KTEB). Early sequence-based phylogeographic studies in babblers [Pomatostomus (58)] reported deep divergence across the CB relative to divergence within the eastern and western regions of the continent. Subsequent multilocus analyses of several avian systems revealed mostly Pleistocene divergences across the CB for congeners (59 ⇓ –61), as well as within species (62, 63). Some studies examining divergence across the region have discovered clines (64) or complex reticulate patterns in the form of introgression for example, in butcherbirds (Cracticus), populations east of the KTEB are introgressed with mtDNA from populations of arid-adapted species to the south that expanded during the last glacial maximum, whereas populations west of the KTEB are not introgressed (65).

Gene tree heterogeneity in multilocus phylogeographic datasets of birds (Red-Backed Fairywren, Malurus melanocephalus Poephila grassfinches Climacteris treecreepers), skinks (Two-Spined Rainbow Skink, C. amax Shaded-Litter Rainbow Skink, C. munda), and mammals (Petrogale rock-wallabies) across northern Australia. (UNE) Map of northern Australia showing the KTEB and CB that separate the KIM, TE, and CY faunas. (B) Cloudograms illustrate topological and branch length variation of gene trees. Violin plots represent the distribution of pairwise sequence divergences across the CB, and black dots indicate mean pairwise sequence divergence, or xy. Red dots and lines are estimates and 95% confidence intervals of population divergence across the KTEB, whereas green dots and lines are estimates and 95% confidence intervals of population divergence across the CB. (C) Distribution of rooted triplets shows that gene trees exhibiting deeper divergence times across the CB than the KTEB are the most frequent in all taxa except the Shaded-Litter Rainbow Skink. Additional details are provided in Texte SI.

Fewer multilocus phylogeographic studies have been conducted for mammals across the monsoonal tropics, yet some common themes emerge. Rock-wallabies (Petrogale), specialists of rocky habitats as the name implies, are strongly structured across the disjunct sandstone plateaus of the region, with deeper divergences across the CB than across the KTEB (66, 67) (Fig. 3). Other macropod species have different degrees of geographic and genetic discontinuity across the CB, suggesting the species’ ecology has played a key role in their ability to adapt and persist across this region. The antilopine wallaroo (Macropus antilopinus), a savanna-woodland specialist, has a disjunct distribution but shallow divergence on either side of the CB grasslands, suggesting recent gene flow or range expansion. By comparison, a more ecologically generalized and widespread congener, the common wallaroo (Macropus robustus), has a more continuous distribution but substantial divergence across the CB (68). Preliminary analyses of other mammals suggest deep divergences across the CB, but require more expansive molecular study (69, 70).

Low-dispersal species, such as lizards and frogs, have been the subject of a burst of recent multilocus phylogeographic studies across the region, including early applications of NGS in this context. Relative to birds and mammals, these taxa exhibit phylogeographic structure at a finer spatial scale, often with cryptic species and greater phylogenetic depth among regions, possibly reflecting a combination of lower dispersal and higher localized persistence through cycles of harsh climate. Deep structure across the CB, and often also the KTEB, is observed across phylogenetically and ecologically diverse reptiles, including species complexes of agamid lizards (71), rainbow skinks (12, 72), several species complexes of geckos (73, 74), and toadlet frogs (75). In many cases, the divergence across the CB appears at deep phylogenetic scales rather than within species. Par exemple, Carlia rainbow skinks have radiated across the KIM and TE, yet these taxa diverged from the eastern species of Carlia in the mid-Miocene (72). Analyses of ∼2,000 exons for the Two-Spined Rainbow Skink (Carlia amax) inferred recent population expansion from western KIM across the KTEB to the western TE, emphasizing that the KTEB is a more porous filter than the CB (12). Studies of low-dispersal taxa (12, 73, 75) are also revealing congruent patterns at a finer scale than the major barriers envisioned by Cracraft (56). These congruent patterns include deep structuring between offshore islands and mainland populations and an unexpected north-south split from the TE to the northern desert region (12, 73, 75). Closely related northern desert taxa often have ranges that are more widespread than those ranges across the savannas and sandstones to the north, and sometimes with evidence of broad-scale introgression (12, 73, 75). For example, in contrast to strong phylogeographic structuring within C. amax, an arid-adapted congener, Carlia munda (Shaded-Litter Rainbow Skink) includes a single widespread clade from the west coast across the northern desert to the east coast (Fig. 3).

Gene Tree Heterogeneity Across the CB.

The complex landscapes and dynamic climate history across this region have resulted in a combination of often strongly vicariant processes across the CB and a mix of divergence and dispersal or introgression across the KTEB. Given that gene tree heterogeneity arises from both ILS and gene flow between populations (Fig. 2), we can expect to see a more dominant phylogenetic signal across the CB in which the deepest split for a majority of gene trees spans the CB, with fewer loci having their deepest split across the KTEB, or between the CY/EF and TE (Fig. 3). We explored this hypothesis for exemplar avian, mammal, and lizard taxa for which we had multilocus sequence data spanning these geographic regions (Fig. 3 and Tables S1 and S2). As expected, among four-tip gene trees (one allele sampled for the KIM, TE, and CY/EF plus outgroup), we found diverse gene tree distributions across the region, with gene trees exhibiting deeper divergence times across the CB than the KTEB being the most frequent (Fig. 3 and Table S3). An exception to this pattern is C. munda, the more arid-adapted lizard, in which the dominant gene tree is one in which the TE and CY alleles are sisters, implying a more isolation-by-distance than vicariance model (76). Analyzing the larger datasets in which these simple gene trees are embedded with coalescent models (77) uniformly suggests deeper population divergence and speciation across the CB than across the KTEB, although these divergences are quite close temporally in several cases (Fig. 3 and Fig. S1). Although our sample sizes are small, the analysis also suggests that the highest genetic diversity currently segregating within each complex varies among regions in Fairy Wrens and wallabies, the highest diversity is in the CY/EF, whereas diversity is similar among regions in C. amax. Finally, the analysis does not support a key prediction of the simplest vicariance scenario: that the effective population sizes of descendant lineages are smaller than the sizes of the ancestral populations inhabiting the area before the vicariant event. Our estimates of ancestral Ne for at least four of the six species groups are smaller than for contemporary lineages in the KIM, TE, or CY. A challenge with our brief analysis of comparative phylogeography across the monsoon tropics of Australia is the diversity of markers, which prevents easy comparison across groups due to differences in substitution rate. Use of a common set of markers, such as provided by various forms of target capture (78 ⇓ –80), whether coding or noncoding, will be an important focus of future research.

Reticulation Driven by Ecology and Introgression.

Our comparative phylogeographic analysis for northern Australia highlights the complex mix of divergence and reticulation and diverse spatial and temporal scales of phylogeographic structure that can emerge. Much of this heterogeneity appears to relate to differences among species in their capacity to persist or disperse across the landscape as climates oscillated over the Quaternary. Whereas it is convenient to focus on common patterns of divergence, in a classic vicariance mindset, closer attention to differing outcomes of reticulation, such as we have seen when comparing the results for C. munda with the results for other taxa spanning the CB, will yield more insight into speciation processes (81).

Following secondary contact, genetically distinct populations can form “tension zones,” maintained over time by a balance between dispersal and selection against hybrids (82), progressively merge via introgression [i.e., ephemeral taxa (46)], or overlap while maintaining their integrity (Fig. 4). A special case of introgression occurs when an expanding lineage overrides a static (relictual) one, but is itself invaded by genes from the resident population due to sequential founder events during the spatial expansion (83). Over time, introgressed chromosome segments will recombine between lineages, leading to a mosaic of coalescent histories within and across loci. These reticulation events can manifest at two scales: in genetic clines, for single-nucleotide polymorphisms (SNPs) at the contact zone(s) themselves, and in lineage-scale migration, as estimated from allopatric populations using isolation-migration (IM) models (Fig. 4). Genome-scale data are enabling new approaches (reviewed in 84), including genomic clines (85) and analyses of lengths of introgressed haplotype blocks (86). Given estimates of recombination rate, the length of immigrant haplotypes can, in principle, be used to estimate the timing of recent introgression events at a lineage scale, a parameter that has proved difficult to infer from IM models (87).

Contrasting processes and views of introgression. (UNE) Progression over time from population splitting, divergence in isolation, and secondary contact, with alternate outcomes: (je) tension zone, (ii) merging, and (iii) overriding of expanding population (blue) over the resident population with introgression from yellow→blue for some genes. (B et C) Contrasting perspectives on introgression among cryptic lineages of Australian Wet Tropics lizards at the local scale (B, contact zone) vs. lineage-scale estimates from IM analyses (modified from ref. 53). Note decreasing introgression at contact zones with increasing divergence time of lineage pairs, but no corresponding signal of decreasing migration at the lineage scale.

In the context of comparative phylogeography, insights into reticulation processes can be gleaned by comparing outcomes for taxa with varying ecologies and lineage divergence times across a common geographic and paleoenvironmental setting. Suture zones can be useful for this purpose, where multiple taxa have co-occurring contact zones (88). The fauna endemic to the rainforest of northeastern Australia are a case in point. Climate-driven fluctuations of rainforest-dependent taxa on mountain tops have resulted in spatially concentrated contact zones between morphologically indistinguishable but genetically distinct lineages (52). A comparative analysis of clines and genetic disequilibria across different contact zones (53) revealed less introgression and stronger genetic disequilibrium between more divergent lineage pairs, showing that reproductive isolation between these phenotypically cryptic lineages scales with divergence time (Fig. 4). However, at the lineage scale, levels of gene flow inferred from IM analyses of comparative transcriptomes are generally low and do not scale with divergence time (Fig. 4). These contrasting patterns remind us that estimates of gene flow are often averaged over the entire divergence history.

The Nexus of Comparative Phylogeography and Speciation Genomics.

So how will a fully genomic perspective enrich our understanding of the nexus between comparative phylogeography and speciation? A plethora of recent whole-genome comparisons among sister taxa reveal fascinating, but complex, heterogeneity of divergence across the genome (reviewed in ref. 84). The most common outcome among recently diverged taxa is stronger differentiation on X and Z sex chromosomes than autosomes and scattered “islands” of high divergence against a background of low divergence. Islands of divergence were initially taken as suggesting locations of incompatible genes in the context of ongoing gene flow (89). However, it is also possible that they reflect varying levels of background selection in the absence of gene flow (90, 91), leading to reinterpretation of some high-profile examples (92).

A key factor emerging from these studies and earlier scans of intraspecific diversity is the strong effect of recombination rate variation on the spatial patterning of genomic diversity, mediated most strongly by hitchhiking (93). Thus, we expect to see reduced within-lineage diversity in regions of low recombination, with a corresponding increase in divergence using measures that are sensitive to levels of within-lineage diversity [e.g., Wright's fixation index Fst (94)]. Paradoxically, it has also been proposed that low-recombination regions, as might occur within chromosomal inversions or near centromeres, will accumulate locally adapted alleles, thereby contributing to genetic incompatibility between lineages (95). Empirical evidence for this proposal is mixed, but there are some positive examples (96 ⇓ –98). Finally, genome comparisons among closely related taxa have also highlighted introgression of adaptive alleles from one lineage to another (35, 99), an old concept reborn (100). Such alleles can readily flow across contact zones even if there is strong hybrid breakdown. Analytical challenges aside, such cases point to the exciting prospect of understanding how adaptive evolution influences divergence and reticulation among lineages.

One limitation of many of the above analyses of genomic divergence during speciation is that the historical biogeographic and environmental setting of isolation and reconnection of diverging lineages over time is often not well established (84). Understanding these processes is the core business of phylogeography, and closer interaction between analyses of historical biogeography and speciation genomics can be expected to bear fruit. Conversely, in the genomic era, comparative phylogeographers will not just have to master details of environmental history, species’ ecology, and the plethora of methods for NGS and demographic inference but will also have to comprehend effects of selection and recombination rate variation across the genome. This challenge is exciting and will serve to strengthen further the link between population genomics and phylogenetics.


Is there a measure for degree of reticulation in a phylogenetic network? - La biologie

When working with phylogenetic networks, how does one talk about reticulate (i.e., web-like) the branches are? Is there a standard measure of this? If so, can someone recommend an R package that can calculate it?

Can you clarify what you mean by "how web-like" a tree is?

A tree can't be web-like since any two nodes are connected by only one path. Maybe just counting the number of internal nodes is what you're after. What is the question you're trying to address with this ?

Not phylogenetic trees, phylogenetic networks. By "web-like" I mean the horizontal lines between branches.

Sorry, I was thrown off by the previous comment. I am not familiar with phylogenetic networks. I would think that various concepts from graph theory could be useful but which to choose would depend on what you're trying to achieve.

Login before adding your answer.

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.


Méthodes

We first introduce the basic concepts and notation, then recap the decomposition technique for arbitrary phylogenetic networks, and finally describe the algorithms for the CCP and the SRF distance.

Concepts and notation

Laisser X be a set of taxa. A rooted phylogenetic network (network for short) over X is an acyclic digraph in which the leaves (i.e., nodes of out-degree zero) are bijectively mapped to X. A taxon typically represents some extant organism or species. A network has a unique root (of in-degree zero).

There can be two types of internal nodes in a network: tree nodes, which include the root and nodes of in-degree one and out-degree of at least one, and reticulation nodes, which have out-degree one and in-degree of at least two. The tree nodes represent speciation events and the reticulation nodes represent reticulation events. We allow degree-two nodes in a network.

Here, we use the following notation for a network N:

T(N): the set of tree nodes in N.

L(N): the set of leaves in N.

R(N): the set of reticulation nodes in N.

V(N): the set of all nodes in N, à savoir T(N) ∪ L(N) ∪ R(N).

E(N): the set of edges in N.

??(N): the root of N.

NE: the subnetwork (V(N),E(N) ∖ E) for a subset EE(N).

NS: the subnetwork (V(N) ∖ V(S),E ′ ), where E ′ =<(X,oui) ∈ E(N) | <X,oui> ⊆ V(N) ∖ V(S)> for a subnetwork S de N.

Pour vous,vV(N), vous est un parent de v et v est un enfant de vous if (vous,v) ∈ E(N). Nous utilisons c(r) to denote the unique child of rR(N). If there is a direct path from vous à v, v s'appelle un descendant de vous.

We use [ r] N to denote the subnetwork below rV(N), which consists of all the descendants of r and the edges between them in N. For a leaf ?? au dessous de r, we use N−[ r] N+?? to denote the subnetwork obtained by replacing [ r] N avec ?? pour que ?? becomes the child of r.

If each reticulation node in a network has exactly two parents, the network is bi-combining. A bi-combining network is binaire if each tree node is of out-degree two. A phylogenetic tree is a binary network without reticulation nodes. If the unique child of each reticulation node in a network is a tree node or a leaf, this network is called réduit.

Following Gunawan et al. [16], we allow a network to have dummy nodes (i.e., unlabelled nodes of out-degree zero) because such a network may be generated in a recursive step of our algorithms.

Given the set of taxa X, une grappe is any proper subset of X (excluding the empty set and the full set). A cluster is trivial if it contains only one element.

In a phylogenetic tree T plus de X, each non-root node induces a unique set of taxa that consists of the labels of the leaves below the node, which is called the cluster of the node. A cluster is in T if it is the cluster of some node in T.

Given a network N plus de X and a phylogenetic tree T plus de X, we say that T est displayed dans N si T can be obtained from N by the following operations: removing all but one incoming edges for each reticulation node in N, removing nodes that are not in any path from ??(N) to a leaf ??X, and contracting degree-two nodes (i.e., nodes of in-degree one and out-degree one). To contract a degree-two node w which has two incident edges (vous,w) et (w,v), we merge the two edges into one edge (vous,v).

A cluster BX est un soft cluster dans N if there is a tree T displayed in N tel que B is a cluster in T. A tree node in a network may represent multiple soft clusters, which could be obtained from different trees displayed in the network. Nous utilisons S C(N) to denote the set of soft clusters in N.

Étant donné BX and a network T au X, the CCP asks whether B is a soft cluster in N [10], which is formulated as below:

Exemple: A phylogenetic network N over a set of taxa X et BX.

Question: Is BS C(N)?

Laisser N 1 et N 2 be two networks over the same set of taxa X. The SRF distance between them is defined to be (|S C(N 1) ∖ S C(N 2)|+|S C(N 2)S C(N 1)|)/2 denoted by SRF(N 1,N 2).

It is worth noting that the SRF distance is not a strict metric, since two distinct networks may represent the same set of soft clusters and hence the SRF distance between them will be zero [10]. Nevertheless, the SRF distance provides a useful measure of network dissimilarity.

A decomposition theorem

The key to solving the CCP and computing the SRF distance is the decomposition theorem, which was first proposed by Gunawan et al. [16] for reticulation-visible networks and used later for arbitrary networks in [20].

The decomposition theorem says that an arbitrary network N can be decomposed into a set of connected tree components which are separated by reticulation nodes. Specifically, there is a tree component C r for each rR(N) ∪ <??(N)>, which is either <c(r)> if rR(N) et c(r) ∈ R(N), or a subtree induced by all the tree nodes and leaves that are reachable from r. A tree component is trivial if it contains only one leaf or if it is empty (for a dummy reticulation node).

A node is visible on a leaf ?? if it lies on all the paths from ??(N) à ??. If a node rR(N) ∪ <??(N)> is visible on a leaf ??, its tree component C r est visible au ?? également. Given two tree components (C_phantom !>) and (phantom !>C_>) , r ′ and (phantom !>C_>) are right below (phantom !>C_) if a parent of r ′ is in (phantom !>C_>) . A tree component is exposé if it contains only one leaf or if all the tree components right below it are trivial.

Évidemment, N contains at least one exposed non-trivial tree component. In addition, an exposed tree component C r is visible if and only if C r contains a leaf or if a reticulation node r ′ exists right below C r such that all the parents of r ′ are in C r.

These concepts are briefly illustrated in Fig. 1. See [16, 20] for more details of the decomposition theorem.

A network N and its tree components. There are nine tree components in N. Five of these components are non-trivial: C r, C r1, C r2, C r5, et C r6, où C r6=. C r7 and r7 are right below C r5, C r2, et C r. C r is visible on all the leaves. C r1 et C r2 are visible, but neither of them is exposed. C r5 is exposed but not visible

Description of the algorithm

The CCP algorithm

With the aid of the generalized decomposition theorem, we extend the linear-time CCP algorithm for reticulation-visible networks in [16] to arbitrary networks.

Roughly speaking, our new CCP algorithm works as follows:

To determine whether or not a cluster C is in a phylogenetic network N, the algorithm selects a non-trivial exposed component M of N. If M is visible, we either find the negative answer to the problem by working on M or we obtain an instance of the problem that is simpler than the input instance (C,N) in linear time proportional to the size of M. In the latter, we reduce the original instance of the CCP to a simpler instance.

If M is not visible, there is then a reticulation node which has a unique leaf child and does not have all parents in M. In this case, two phylogenetic networks N 1 et N 2 are derived from N, which contain fewer nodes than N. The algorithm is then called on both instances ( C,N 1 ) et ( C,N 2 ) recursively.

Although this algorithm seems simple, it has significantly less time complexity when the input network is binary. In the rest of this section, we present a formal description of the algorithm.

Laisser N be a network over X et BX, respectivement. We examine a non-trivial exposed tree component C r de N.

The reticulation nodes below C r are divided into inner-reticulation nodes for which the parents are all in C r, and cross-reticulation nodes for which some parents are not in C r. Nous utilisons je R(C r) et C R(C r) to denote the sets of inner- and cross- reticulation nodes, respectively. For example, in Fig. 1, je R(C r5)= et C R(C r5)=.

Nous utilisons L r to denote the set of leaves on which C r is visible:

We use (check _) to denote the set of leaves below C r which are in B and on which C r is not visible:

For example, in Fig. 1, L r5= and we can get (check _< ext >=< ext , exte >) when assuming B=.

Supposer que L r is non-empty. C r is then visible with respect to a leaf ??L r. We first check whether B is a soft cluster in C r. This can be solved by a linear-time algorithm [16]. If not, we then solve the CCP according to the relationship between L r et B.

Let (ar =X ackslash B) . Si L rB and ( L_ cap ar eq emptyset ) , B must be a soft cluster of a node in C r si B is a soft cluster in N [16].

If (L_ cap ar = emptyset ) , B may be a soft cluster of ??(C r) or a node in a larger subnetwork containing C r. Assuming that r ′ ∈ C R(C r), we then define:

Depuis L rB and (check _ subseteq B) , (hat subseteq B) . If (hat = B) , B is a soft cluster of ??(C r) dans N une. Otherwise, if (hat subset B) , we set:

Si L rB= , B may be a soft cluster of a node in C r if (check _ eq emptyset ) . Otherwise, when B is not a soft cluster of a node in C r et r ′ ∈ C R(C r), we define:

With this notation, we can get Theorem 1 for arbitrary networks, which is similar to a theorem proved for reticulation-visible networks in [16]. Theorem 1 is proved in the Additional file 1.

Theorem 1

Suppose que C r is a non-trivial, exposed and visible tree component in a network N over the taxa set X, et cela BX. Laisser L r, (hat ) , B ′ , N une′, and N b′ be defined above.

(je) If (hat subset B) , B is a soft cluster in N si et seulement si B ′ is a soft cluster in N une′.

(ii) Si B is not a soft cluster of a node in C r et L rB= , B is a soft cluster in N si et seulement si B is a soft cluster in N b′.

Supposer que C r is not visible. Si C r≠<c(r)>, there is at least one reticulation node r ′ right below C r such that (phantom !>C_) is trivial and at least one parent of r ′ is not in C r. Si C r=<c(r)> and (c(r)=r'phantom !>) , then at least one parent of r ′ is not r. We can now define: