Informatie

Is er een maat voor de mate van verknoping in een fylogenetisch netwerk?


Als je met fylogenetische netwerken werkt, hoe praat je dan over netvormige (d.w.z. webachtige) takken? Is hier een standaard maatstaf voor?


Genealogieën: stambomen en fylogenieën zijn netwerken van netwerken, niet alleen uiteenlopende bomen

Stambomen illustreren de genealogische relaties tussen individuen, en fylogenieën doen hetzelfde voor groepen organismen (zoals soorten, geslachten, enz.). Hier geef ik een kort overzicht van de huidige concepten en methoden voor het berekenen en weergeven van genealogische relaties. Het is al lang bekend dat deze relaties verknoopt zijn in plaats van strikt uiteenlopend, en dus worden zowel stambomen als fylogenieën correct behandeld als netwerken in plaats van bomen. Momenteel worden de meeste stambomen in plaats daarvan gepresenteerd als "stambomen", en de meeste fylogenieën worden gepresenteerd als fylogenetische bomen. Desalniettemin toont de historische ontwikkeling van concepten aan dat netwerken ouder waren dan bomen op de meeste gebieden van de biologie, inclusief de studie van stambomen, biologietheorie en biologiepraktijk, evenals in historische taalkunde in de sociale wetenschappen. Bomen zijn eigenlijk geïntroduceerd om een ​​eenvoudiger conceptueel model voor historische relaties te bieden, aangezien bomen een specifiek type eenvoudig netwerk zijn. Computationeel maken bomen en netwerken deel uit van de grafentheorie, bestaande uit knooppunten die zijn verbonden door randen. In deze wiskundige context verschillen ze alleen in de afwezigheid of aanwezigheid van reticulatieknooppunten, respectievelijk. Er zijn twee soorten grafieken die fylogenetische netwerken kunnen worden genoemd: (1) gewortelde evolutionaire netwerken en (2) onbewortelde affiniteitsnetwerken. Er zijn nogal wat rekenmethoden voor niet-gewortelde netwerken, die twee hoofdrollen hebben in de fylogenetica: (a) ze fungeren als een generieke vorm van multivariate gegevensweergave en (b) ze worden specifiek gebruikt om haplotype-netwerken weer te geven. Evolutionaire netwerken zijn moeilijker af te leiden en te analyseren, omdat er geen wiskundig algoritme is om unieke historische gebeurtenissen te reconstrueren. Er is dus momenteel geen coherent analytisch kader voor het berekenen van dergelijke netwerken.

Dit is een voorbeeld van abonnementsinhoud, toegang via uw instelling.


Voorrondes

Gegeven een verzameling taxa X . Een subset van X (behalve zowel ∅ als X ) wordt een cluster genoemd. Voor twee clusters C1 en C2 op X , als ze onsamenhangend zijn of als de ene een subset is van de andere, d.w.z. C1C2 = ∅ of C1C2 of C2C1, we zeggen dat C1 en C2 compatibel zijn, anders zijn ze incompatibel. Gegeven een verzameling clusters C op X . Als twee clusters in C compatibel zijn, wordt C compatibel genoemd, anders is het incompatibel. De incompatibiliteitsgrafiek IG( C ) = (V, E) van C wordt gedefinieerd als een ongerichte graaf met knooppuntenset V = C en randset E, waarbij twee clusters zijn verbonden door een rand als en alleen als ze incompatibel zijn.

Laten S een deelverzameling van X zijn. De beperking van C tot S, aangeduid met C | S , wordt gedefinieerd als het resultaat van het verwijderen van alle taxa in X S (d.w.z. de taxa in X maar niet in S) van elk cluster in C . De incompatibiliteitsgraad van C, aangegeven met d(C), wordt gedefinieerd als het aantal randen in I G(C). Laat C een verzameling clusters zijn op X en x een taxon in X . De incompatibiliteitsgraad van x, aangeduid met NS(x), wordt gedefinieerd als het resultaat van het aftrekken van de incompatibiliteitsgraad van C | X < x >van die van C , d.w.z. d ( x ) = d ( C ) - d ( C | X < x >) . Als de incompatibiliteitsgraad van een taxon x is maximaal onder alle taxa in X , d.w.z. d ( x ) = max < d ( y ) | y ∈ X >, dat noemen we x is het incompatibiliteitstaxon w.r.t. C . De frequentie van een taxon x w.r.t. C , aangeduid met F(x), wordt gedefinieerd als het aantal clusters met x, d.w.z. f ( x ) = | < C ∈ C | x ∈ C >| .

Een geworteld fylogenetisch netwerk N = (V, E) op X is een gewortelde gerichte acyclische graaf (kortweg DAG), en de bladeren zijn bijectief gelabeld als X . De ingraad van een knoop vV wordt aangeduid als σ(v). een knoop v is een netvormige knoop als σ(v) ≥ 2 anders is het vooral een boomknooppunt, een boomknooppunt is een wortelknooppunt als σ(v) = 0. Een rand e = (jij, v) is een boomrand als v is een boomknooppunt, anders is het een netvormige rand. Het reticulatienummer in een netwerk N = (V, E) is ∑σ(v(σ(v) − 1) = |E| − |V| + 1.

Gegeven een cluster C en een gewortelde fylogenetische boom t. Als er een rand is e in t zodat de verzameling van taxa bereikbaar is vanaf e gelijk aan C, we zeggen dat t vertegenwoordigt C. Gegeven een netwerk N, bij het verbinden van één inkomende rand en het loskoppelen van alle andere inkomende randen voor elk netvormig knooppunt, als er een boomrand bestaat e zodat de verzameling van taxa bereikbaar is vanaf e gelijk aan C, we zeggen dat N vertegenwoordigt C in de zachte zin. Als alternatief, als er een boomrand is e in N zodat de verzameling van taxa bereikbaar is vanaf e gelijk aan C, we zeggen dat N vertegenwoordigt C in de bekabelde zin.

Laten N = (V, E) een netwerk zijn dat de verzameling clusters C vertegenwoordigt. Een cluster C ∈ C wordt vaak vertegenwoordigd door meer dan één boomrand in N en een boomrand eE vertegenwoordigt vaak meer dan één cluster in C . Als er een toewijzing bestaat ϵ van C naar de verzameling boomranden van N, waar ϵ(C) is een boomrand die staat voor C voor C ∈ C , zodat voor elke twee clusters C 1 , C 2 ∈ C , C1 en C2 in dezelfde verbonden component van de incompatibiliteitsgrafiek liggen IG( C ) als en slechts als ϵ(C1) en ϵ(C2) zijn opgenomen in dezelfde biconnected component van N. Dan noemen we dat N is afbreekbaar tov C .

Bij het construeren van geroote fylogenetische netwerken van geroote fylogenetische bomen, berekenen de methoden eerst de clusters die worden vertegenwoordigd door de invoerbomen en construeren vervolgens een geroot fylogenetisch netwerk dat alle clusters vertegenwoordigt.

De gewortelde fylogenetische netwerken kunnen de evolutionaire geschiedenis beschrijven in de aanwezigheid van netvormige gebeurtenissen, zoals horizontale genoverdrachten, hybridisaties en recombinaties. Deze netvormige gebeurtenissen zijn in werkelijkheid zeldzaam [33]. Dienovereenkomstig wordt verwacht dat het geconstrueerde netwerk het minimale aantal netvormige knooppunten heeft. Laten N een netwerk zijn dat is geconstrueerd voor de invoerclusterset C . Neem aan dat C′ de verzameling clusters is die wordt vertegenwoordigd door N. In feite bevat C′ meer clusters dan C , d.w.z. C ⊊ C ′ . Hier definiëren we de redundante clusters C 0 van N als de clusters die wel in C ′ zitten maar niet in C , d.w.z. C 0 = C ′ C . In fylogenetische analyse zijn de taxa in een cluster vermoedelijk monofyletisch. Bijgevolg zou de ideale situatie Co = ∅ zijn, d.w.z. alle clusters die door het geconstrueerde netwerk worden vertegenwoordigd, zouden de clusters zijn die in de invoerbomen worden weergegeven, en geen andere. Daarom, door middel van het spaarzaamheidsprincipe, is het best geconstrueerde netwerk er een dat het aantal redundante clusters minimaliseert, wat gebaseerd is op de voorwaarde dat het een minimaal aantal netvormige knooppunten heeft.

De volgende sectie introduceert de belangrijkste methoden voor het construeren van gewortelde fylogenetische netwerken van gewortelde fylogenetische bomen.


Resultaten

Waarschijnlijkheidsdichtheid van een genboom.

Gegeven een fylogenetisch netwerk Ψ en een genealogie G j voor locus J (topologie en vertakkingslengtes in beide gevallen), duiden we met H Ψ ( G j ) de verzameling van alle coalescentiegeschiedenissen van G j binnen de takken van Ψ aan. Vervolgens wordt de verdeling (dichtheid) van genenbomen gegeven door p ( G j | Ψ , Γ ) = ∑ h ∈ H Ψ ( G j ) p ( h | Ψ , Γ ) , [4] waarbij Γ de overervingskansen is matrix, zoals hierboven beschreven. Voor een rand b = ( x , y ) ∈ E ( ) definiëren we T b ( h ) als de vector van tijden (in oplopende volgorde) van samensmeltingsgebeurtenissen die plaatsvinden op vertakkingen B onder de samenvloeiende geschiedenis H en de tijd van de knoop ja (een formele definitie wordt gegeven in SI-bijlage). We geven aan met T b ( h ) [ i ] de lhet element van de vector. Verder geven we met u b ( h ) het aantal genlijnen aan dat edge binnenkomt B en door v b ( h ) het aantal genlijnen dat edge verlaat B onder H. Dan hebben we p ( h | Ψ , Γ ) = ∏ b ∈ E ( Ψ ) [ ∏ i = 1 | Tb ( h ) | − 1 e ( ub ( h ) − ik + 1 2 ) ( T b ( h ) ik + 1 - T b ( h ) ik ) ] × e ( vb ( h ) 2 ) ( τ Ψ ( b ) T b ( h ) | T b ( h ) | ) × Γ [ b , j ] ub ( h ) , [5] waarbij τ Ψ ( b ) voor rand b = ( x , y ) is de tijd van de knoop x in het fylogenetische netwerk Ψ. Een volledige afleiding van de formule en een efficiënter algoritme om het te berekenen in de trant van Yu et al. (17), die expliciete sommaties over de mogelijke samengevoegde geschiedenissen vermijden, worden gegeven in: SI-bijlage.

Zoeken in de ruimte van fylogenetische netwerken.

Laten we Ω ( n ) de ruimte aanduiden van alle fylogenetische netwerken op N taxa duiden we met Ω ( n , k ) de deelruimte van Ω ( n ) aan die alle fylogenetische netwerken (rDAG's) bevat met N bladeren en k reticulatie knooppunten. In het bijzonder is Ω ( n , 0 ) de deelruimte die alle fylogenetische bomen bevat. Om de fylogenetische netwerkruimte op een gelaagde manier te doorzoeken, definiëren we twee operaties die het mogelijk maken om binnen Ω ( n , k ) te zoeken naar een gegeven k: een bewerking waarmee het zoeken een laag van Ω ( n , k ) naar Ω ( n , k + 1 ) kan laten stijgen en een bewerking waarmee het zoeken een laag van Ω ( n , k ) naar Ω ( n , k 1 ). Voor het zoeken binnen een laag, verplaatsen de bewerkingen ofwel de bestemming van een vernettingsrand of verplaatsen de bron van een rand (vernetting of niet). Voor het opstijgen van een laag bestaat de bewerking uit het toevoegen van een vernettingsrand tussen twee bestaande randen in het netwerk, en voor het afdalen van een laag verwijdert de bewerking een vernettingsrand (meer details vindt u in SI-bijlage). Het is vermeldenswaard dat, hoewel de ruimte van alle fylogenetische boomtopologieën op N taxa is eindig, de ruimte van alle fylogenetische netwerktopologieën op N taxa is in theorie oneindig, omdat Ω ( n ) = ∪ k ≥ 0 Ω ( n , k ) en k zijn onbegrensd. Beschouw bijvoorbeeld het geval van slechts twee taxa. Er is in dit geval een unieke, gewortelde boom. Omdat er echter meerdere hybridisatiegebeurtenissen kunnen plaatsvinden tussen dezelfde twee zustertaxa op verschillende tijdstippen, kan een willekeurig aantal horizontale randen tussen deze twee taxa worden toegevoegd. Desalniettemin, of dergelijke repetitieve hybridisatiescenario's identificeerbaar zijn uit typische genomische datasets, is een andere vraag.

Een heuristiek voor het schatten van de optimale vertakkingslengtes voor een boomtopologie met een vaste soort, gegeven genboomtopologieën, die is gebaseerd op herhaalde toepassing van de methode van Brent (22) werd geïntroduceerd door Wu (20). We gebruiken een vergelijkbare heuristiek voor het schatten van de fylogenetische netwerktaklengtes en overervingskansen (volledige details worden gegeven in SI-bijlage). Door topologische transformaties en parameterschattingsheuristieken te koppelen aan de waarschijnlijkheidsformulering hierboven, kan de ruimte op een heuvelklimmende manier worden doorzocht om een ​​ML-fylogenetisch netwerk af te leiden. Gezien het bestaan ​​van lokale optima binnen elke laag, kunnen meerdere, onafhankelijke runs worden gemaakt.

Controleren op modelcomplexiteit.

Omdat netwerken in Ω ( n , k + 1 ) complexere modellen bieden dan netwerken in Ω ( n , k ) , moeten de hierboven beschreven benaderingen het modelselectieprobleem aanpakken. Informatiecriteria zijn al gebruikt in de context van fylogenetische netwerken (12, 16), en we gebruiken ze hier (in plaats van te zoeken op basis van de waarschijnlijkheidsscore, wordt er gezocht op basis van de waarden van deze criteria, waarin de waarschijnlijkheidsscores zijn verwerkt) . Een andere benadering die we hier voorstellen, voor het eerst voor zover wij weten, is het gebruik van K-fold cross-validatie, waarbij de input set van genenbomen wordt opgedeeld in K subsets van gelijke grootte, worden de parameters van het model afgeleid uit K 1 subsets en wordt de pasvorm van het model van de resterende subset berekend. Deze fit wordt berekend door de frequenties van de genenbomen in de validatiesubset te vergelijken met de verdeling van de genbomen die door het afgeleide netwerk worden geproduceerd. Als de fit van het beste netwerk Ψ ″ gevonden in Ω ( n , k + 1 ) niet veel beter is (we hanteren een afkapwaarde van 3% verbetering, gekozen op basis van empirische waarnemingen) dan de fit van het beste netwerk Ψ ′ gevonden in Ω ( n , k ), we verklaren k om de juiste schatting te zijn van het aantal verknopingsknopen en Ψ ′ om het optimale fylogenetische netwerk te zijn. Het is belangrijk op te merken dat dit idee van kruisvalidatie alleen werkt voor volledig opgeloste genboomtopologieën, omdat in het geval van genenbomen met vertakkingslengtes de frequenties van de genbomen in de validatiesubset niet informatief zijn.

Ten slotte, om de ondersteuning van de fylogenetische netwerken die we afleiden te beoordelen, stellen we voor om parametrische bootstrapping te gebruiken. Nadat we een netwerk Ψ hebben afgeleid uit de data G , gebruiken we Ψ om ℓ datasets te genereren, waaruit we ℓ fylogenetische netwerken Ψ 1 , … , Ψ ℓ afleiden. We schatten dan de ondersteuning van elke tak B in Ψ als het aantal netwerken (van de ℓ ) met een equivalente vertakking. We zeggen dat twee randen in twee fylogenetische netwerken equivalent zijn als (l) een of beide zijn vernettingsranden of beide niet en (ii) beide definiëren dezelfde clusters (het cluster gedefinieerd door een vertakking is de verzameling van alle taxa onder die vertakking in het netwerk).

Prestaties op gesimuleerde gegevens.

We hebben alle hierboven beschreven methoden geïmplementeerd in het openbaar beschikbare, open-source softwarepakket PhyloNet (23) en de prestaties van de methoden op verschillende gesimuleerde datasets bestudeerd. In de simulatiestudie waarvan de resultaten worden gerapporteerd in Fig. 2, gebruikten we het fylogenetische netwerk Ψ 1 als het modelnetwerk, en voor verschillende aantallen loci evolueerden we genbomen onder de coalescentie binnen de takken van het netwerk en simuleerden vervolgens sequentie-evolutie op deze genenbomen met variërende sequentielengtes. Vervolgens hebben we voor elke sequentie-uitlijning 100 genenbomen geschat met behulp van ML met bootstrapping. Ten slotte hebben we netwerken afgeleid met behulp van onze ML-methode van (l) echte genboomtopologieën, (ii) geschatte genboomtopologieën, (iii) ware genboomtopologieën en vertakkingslengtes, en (NS) geschatte genboomtopologieën en vertakkingslengtes. De resultaten van (l) en (ii) worden getoond in Fig. 2B, terwijl de resultaten van (iii) en (NS) worden getoond in Fig. 2C. Voor elke instelling van het aantal loci en de lengte van de sequentie hebben we 30 datasets gegenereerd en conclusies getrokken over alle.

Nauwkeurigheid van de methode op gesimuleerde gegevens. (EEN) Gegevens werden gegenereerd langs het fylogenetische netwerk Ψ 1 (alle interne takken, behalve de horizontale rand, hebben een lengte van 1 coalescentie-eenheid en de overervingskans is 0,1 voor alle loci). Resultaten gebaseerd op schattingen van genboomtopologie (B) en gene tree topologie en vertakking lengte schattingen (C) zijn getoond. Voor elk aantal loci komt de meest rechtse balk overeen met gevolgtrekking uit de echte genealogieën en de andere drie balken, van links naar rechts, komen overeen met geschatte genealogieën (met behulp van 100 bootstrap-replica's en Vgl. 3) uit reeksen met een lengte van respectievelijk 250, 500 en 1.000. De donkerblauwe, cyaan en gele gebieden komen overeen met het aantal keren dat elk van de netwerken Ψ 1 , Ψ 2 en Ψ 3 respectievelijk in EEN werd afgeleid. De kastanjebruine regio komt overeen met het aantal keren dat een ander netwerk met een enkele reticulatie werd afgeleid. Hier werd voor elk van de loci één individu bemonsterd per taxon.

Terwijl de hybridisatie in het modelnetwerk B en de meest recente gemeenschappelijke voorouder (MRCA) van C en D omvat, kan de lengte van de vertakking tussen de hybridisatiegebeurtenis en de divergentie van C en D van hun MRCA een groot effect hebben op het onderscheid tussen de echte hybridisatiescenario's en de twee gegeven door Ψ 2 en Ψ 3 in Fig. 2EEN. Daarom hebben we voor elke dataset vastgelegd of de methode een van de drie netwerken afleidde die in figuur 2 worden getoondEEN, in tegenstelling tot elk ander netwerk met een enkele reticulatie.

Verschillende trends kunnen worden waargenomen in Fig. 2EEN. Ten eerste resulteert het gebruik van de echte genboomtopologieën met vertakkingslengten in nauwkeuriger gevolgtrekkingen dan het gebruik van alleen genboomtopologieën. Deze bevinding is niet verwonderlijk, omdat het eerste type gegevens meer informatie bevat dan het laatste. In het bijzonder, wanneer 80 of 160 loci worden gebruikt, is het afgeleide netwerk van de echte genenbomen met vertakkingslengten altijd het echte netwerk. Aan de andere kant, wanneer alleen de genboomtopologieën voor 160 loci werden gebruikt, leverde de gevolgtrekking in vijf van de 30 gevallen een van de twee alternatieve netwerken Ψ 1 en Ψ 2 op. Ten tweede verbeterde de nauwkeurigheid van de gevolgtrekkingen naarmate het aantal loci toeneemt en de sequentielengte toeneemt, hoewel de toename van het aantal loci veel meer een positief effect had op de nauwkeurigheid van de gevolgtrekkingen. Ten derde is een zeer verrassend resultaat dat wanneer alleen genboomtopologieën worden gebruikt, het gebruik van de echte genenbomen bijna nooit resulteert in een betere nauwkeurigheid dan bij het gebruik van geschatte genboomtopologieën voor een bepaald aantal loci. Dit resultaat getuigt van het feit dat de methode, wanneer zorgvuldig rekening wordt gehouden met onzekerheid in de schattingen van de genenboom, zeer goede resultaten kan opleveren. Zelfs bij het gebruik van genboomtopologieën en vertakkingslengtes, is de winst in nauwkeurigheid bij het gebruik van de echte genenbomen erg klein in vergelijking met het gebruik van de genboomschattingen waarbij rekening wordt gehouden met onzekerheid. Ten vierde resulteert de combinatie van een lage waarde van overervingskans (0,1 in deze simulatie) en een relatief korte tijd tussen hybridisatie en daaropvolgende soortvorming in onzekerheid bij het identificeren van de donor en ontvanger van de hybridisatiegebeurtenis. Als u bijvoorbeeld alleen genboomtopologieën gebruikt voor 160 loci, is het afgeleide netwerk altijd een van de drie netwerken Ψ 1 , Ψ 2 en Ψ 3 , zelfs als u denkt dat het meestal Ψ 1 is. We ontdekten dat het vergroten van de vertakkingslengtes of de overervingskansen zou resulteren in hogere nauwkeurigheden. Bovendien vonden we in onze simulaties dat het verhogen van het aantal individuen dat per taxon werd bemonsterd, zou resulteren in een verbeterde nauwkeurigheid, zij het nogal licht (SI-bijlage). We verwachten echter dat het bemonsteren van meer individuen zou resulteren in meer significante verbeteringen op grotere of complexere datasets. In termen van de afgeleide overervingskansen resulteerden de echte genenbomen in zeer nauwkeurige schattingen, terwijl geschatte genbomen met vertakkingslengtes in het algemeen resulteerden in nauwkeurigere schattingen van de kansen. Ten slotte ontdekten we dat kruisvalidatie over het algemeen beter presteert dan informatiecriteria bij het bepalen van het aantal hybridisatiegebeurtenissen (inclusief op de biologische dataset, zoals hieronder besproken). Uitgebreidere simulatieresultaten onder scenario's die gemakkelijker kunnen worden afgeleid dan degene die we hier hebben besproken, zijn opgenomen in: SI-bijlage.

Analyse van een multilocus huismuis dataset.

We hebben onze methode ook gebruikt om een ​​multilocus-dataset van huismuis (M. musculus) genomen, verkregen uit de studies van Staubach et al. (7), Didion et al. (24) en Yang et al. (25) (meer details vindt u in SI-bijlage). Staubach et al. (7) vond substantieel genoom-breed bewijs van subspecifieke introgressie in alle vier de populaties, wat neerkomt op 3% van het genoom in de twee M.m. domesticus populaties (een uit Frankrijk en de andere uit Duitsland), 4% in een M.m. spierpijn bevolking uit Kazachstan, en 18% in an M.m. spierpijn bevolking uit Tsjechië. Het is echter belangrijk op te merken dat de HAPMIX-methode (26), die werd gebruikt door Staubach et al. (7), houdt niet expliciet rekening met ILS.

Onze studie omvatte alle monsters in de studie van Staubach et al. (7). Bovendien omvatte onze studie extra monsters van een M.m. spierpijn populatie uit China (25) die niet zijn gebruikt in de studie van Staubach et al. (7). In deze analyse hebben we alleen geschatte genboomtopologieën gebruikt. De reden hiervoor is dat de genomische sequenties worden verkregen van zeer nauw verwante individuen (deze individuen zijn vijf individuen van dezelfde soort), en er is zeer weinig variatie in de gegevens om de lengte van takken met enige nauwkeurigheid te schatten. Bovendien maakt deze lage variatie een goede bootstrap-analyse van genenbomen voor de individuele loci niet mogelijk. Het krachtige signaal in deze dataset komt van het zeer grote aantal loci. In onze analyse vonden we een significante verbetering in een fylogenetisch netwerk met een enkele reticulatie ten opzichte van geen reticulaties en een significante verbetering in een fylogenetisch netwerk met twee reticulaties ten opzichte van een enkele reticulatie. Toen we echter doorgingen met het zoeken naar het optimale netwerk met drie reticulaties, ontdekten we dat de verbetering die werd verkregen door een derde reticulatiegebeurtenis te beschouwen, niet significant was op basis van de informatiecriteria, en dat er helemaal geen verbetering was op basis van kruisvalidatie. We noemden dus het optimale fylogenetische netwerk met twee reticulaties als onze hypothese voor de evolutionaire geschiedenis van deze reeks genomen. Deze evolutionaire geschiedenis wordt getoond in Fig. 3 (meer details van de resultaten en analyses worden gegeven in SI-bijlage). Het fylogenetische netwerk is niet ultrametrisch en het is de moeite waard om te benadrukken dat de vertakkingslengtes worden gegeven in coalescentie-eenheden. Het gebrek aan ultrametriciteit kan dus te wijten zijn aan verschillende populatiegroottes of, in mindere mate, verschillende generatietijden.

Optimaal fylogenetisch netwerk afgeleid van de huismuis (M. musculus) gegevensbestand. Een enkel individu werd bemonsterd uit elk van de vijf populaties: M.m. domesticus uit Frankrijk (DF), M.m. domesticus uit Duitsland (DG), M.m. spierpijn uit Tsjechië (MZ), M.m. spierpijn uit Kazachstan (MK), en M.m. spierpijn uit China (MC). De analyse vond meerdere, bijna even optimale, fylogenetische netwerken met twee reticulatiegebeurtenissen. Deze meerdere netwerken waren het allemaal eens over de ontvangende populaties, maar waren het niet eens over de donorpopulaties. Eén hybridisatie (de bovenste gestippelde horizontale pijl) omvat de MRCA van DF en DG als een ontvangende populatie, maar lijkt MK, MC of hun MRCA als de donorpopulatie te hebben betrokken. De tweede hybridisatie (de onderste gestippelde horizontale pijl) omvat MZ als een ontvangende populatie, maar lijkt DF, DG of hun MRCA als de donorpopulatie te hebben betrokken. Taklengtes in coalescentie-eenheden (op de boomtakken) en overervingskansen (op de horizontale randen) worden getoond (volledige details van de gegevens en resultaten worden gegeven in SI-bijlage).

Onze analyse van het genoom van huismuizen produceert een evolutionaire geschiedenis die verschilt van die gerapporteerd door Staubach et al. (7) niet alleen in termen van het aantal betrokken populaties, maar ook door rekening te houden met de evolutionaire geschiedenis van de betrokken populaties. We beschouwen de percentages van het genoom met introgressieve oorsprong gerapporteerd door Staubach et al. (7) overschat, omdat introgressie waarbij een voorouderlijke populatie betrokken was die later in meer dan één bestaande populatie werd opgesplitst, in het geval van de studie van Staubach et al. voor elke bestaande populatie meervoudig zou worden gerapporteerd. (7). Aan de andere kant zouden dezelfde percentages worden onderschat in het geval dat vermengde populaties werden gebruikt in plaats van de niet-vermengde referentiepopulaties vereist door HAPMIX, zoals Staubach et al. (7) deed door gebruik te maken van vermoedelijk introgressieve muismonsters om de referentiepopulaties te construeren. Met name vereist onze methodologie niet het gebruik van niet-vermengde referentiepopulaties.

We veronderstellen dat de meer recente introgressie-gebeurtenis in Fig. 3 te wijten is aan de genenstroom van secundair contact, waarbij het bereik van de twee M. musculus ondersoorten overlapten elkaar, ruwweg op de grens tussen Duitsland en Tsjechië. De biologische interpretatie van de oudere introgressiegebeurtenis is minder duidelijk. We vermoeden dat de gebeurtenis verband houdt met genenstroom tijdens en na subspecifieke divergentie. Verder onderzoek kan belangrijke aanwijzingen opleveren voor de mechanistische basis van de evolutie van ondersoorten in M. musculus en het proces van soortvorming zelf.

Het is belangrijk op te merken dat hoewel we een zeer groot aantal loci gebruikten, er nog steeds onzekerheid was in de afgeleide oorsprong van de twee hybridisatiegebeurtenissen (zoals getoond in Fig. 3), een vergelijkbaar patroon als dat waargenomen in de simulatieresultaten en hierboven besproken. Deze onzekerheid is een weerspiegeling van het zwakke signaal in deze gegevens, gekoppeld aan de lage overervingswaarschijnlijkheid en de korte vertakkingslengte tussen de hybridisatie en de MRCA van M.m. spierpijn uit China en M.m. spierpijn uit Kazachstan en M.m. domesticus uit Frankrijk en M.m. domesticus uit Duitsland, wat een probleem is dat we hierboven hebben besproken in de context van de gesimuleerde gegevens. De gebruikte monsters zijn zeer nauw verwant, wat resulteert in genomen met een zeer klein aantal segregerende plaatsen, en dus een zwakker signaal voor gevolgtrekking. Desalniettemin is de onzekerheid gelokaliseerd in die zin dat de potentiële donoren van het genetische materiaal van elke hybridisatie-gebeurtenis rond een enkele voorouderlijke knoop draaien. Omdat alle vijf onderzochte populaties nauw verwant zijn, waren de meeste gereconstrueerde genenbomen niet binair, vanwege identieke sequenties van meerdere allelen. Omdat bootstrapping in deze scenario's niet nuttig is (elke locus heeft een handvol sites, waarvan de meeste monomorf zijn), gebruikten we de niet-binaire genboomtopologieën voor de loci en beschouwden we de verzameling van alle resoluties als de set van te gebruiken genboomschattingen in verg. 3.


Achtergrond

Fylogenetische netwerken zijn grafieken die worden gebruikt voor het weergeven van fylogenetische relaties tussen verschillende taxa, en worden meestal gebruikt wanneer een boomweergave niet voldoende is. Er zijn veel verschillende soorten fylogenetische netwerken en het is nuttig om onderscheid te maken tussen twee hoofdklassen: impliciet fylogenetische netwerken die hulpmiddelen bieden om incompatibele fylogenetische signalen te visualiseren en analyseren, zoals gesplitste netwerken [1, 2], en expliciet fylogenetische netwerken die expliciete scenario's van reticulaire evolutie bieden, zoals hybridisatienetwerken [3-7], HGT-netwerken [8] en recombinatienetwerken [9-19].

De software die momenteel beschikbaar is voor de berekening en analyse van expliciete fylogenetische netwerken bestaat uit een reeks basisimplementaties van algoritmen die zijn ontwikkeld om de rekentaak op te lossen [6, 14, 16-18, 20]. Het grootste deel van de software is commandoregelgestuurd en een aansprekende visualisatie van de resultaten ontbreekt vaak. Het is essentieel om een ​​tool te hebben die zowel een breed gebruik van de methoden waarover biologen ter beschikking staan, als een betere en verdere ontwikkeling van nieuwe methoden mogelijk maakt.

SplitsTree is een applicatie ontwikkeld in onze onderzoeksgroep, oorspronkelijk gericht op de fylogenetische analyse van sets van splitsingen. De nieuwste versie van SplitsTree [2] bevat een verscheidenheid aan methoden voor de berekening, visualisatie en interpretatie van fylogenetische bomen en impliciete fylogenetische netwerken. Twee belangrijke voordelen van SplitsTree zijn de grafische gebruikersinterface (GUI) en de integratie van algoritmen via een interfacegestuurde class loader (plug-ins). In dit artikel presenteren we een uitbreiding van SplitsTree waarmee het programma expliciete fylogenetische netwerken kan verwerken. De uitbreiding lost twee belangrijke problemen op: een efficiënte integratie van expliciete fylogenetische netwerken en het visualiseren van deze netwerken.


Methoden:

Liu et al. heeft onlangs MP-EST geïntroduceerd, een maximale pseudo-waarschijnlijkheidsbenadering voor het schatten van soortbomen uit een verzameling gewortelde genbomen onder de multispecies-coalescentie [26]. De methode resulteerde in significante verbeteringen in de looptijd van statistische inferentie van soortenbomen. Geïnspireerd door dit werk stellen we een methode voor voor het schatten van soortenfylogenieën in de aanwezigheid van zowel hybridisatie als onvolledige afstammingssortering onder maximale pseudo-waarschijnlijkheid.

Fylogenetische netwerken, genenbomen en gewortelde triples

Een (binair) fylogenetisch netwerk [29] Ψ op verzameling X van soorten (taxa) is een gewortelde, gerichte, acyclische graaf waarvan de knooppuntenverzameling V (Ψ) = <R> V L V t V N waar

R is de wortel van Ψ en voldoet aan NS (R) = 0 en NS + (R) = 2

V L: de leaf-set bijectief gelabeld door X , waarbij NS (v) = 1 en NS + (v) = 0 voor willekeurig v V L

V t: interne boomknooppunten, waar NS (v) = 1 en NS + (v) = 2 voor elke v VT en,

V N: verknopingsknooppunten, waar NS (v) = 2 en NS + (v) = 1 voor elke v V N.

Hier, NS (v) en NS + (v) zijn de in-graad en uit-graad van knoop v, respectievelijk. We duiden aan met E(Ψ) de reeks randen in netwerk Ψ. Het fylogenetische netwerk heeft vertakkingslengtes λ : E(Ψ) + . Hierna zullen we Ψ gebruiken om zowel de topologie als de vertakkingslengtes van een fylogenetisch netwerk aan te duiden. Verder, zoals in [22, 24], is er voor een probabilistische instelling een extra functie, de overervingskans genoemd, γ : E(Ψ) [0, 1] die voldoet aan:

γ(e) = 1 voor elke rand e wiens hoofd een boomknooppunt is, en

γ(e1) + γ(e2) = 1 voor elk paar randen e1 en e2 waarvan de kop hetzelfde verknopingsknooppunt is.

In [24] hebben we besproken hoe we de functie kunnen generaliseren γ zodat het per loci varieert, en dat generalisatie triviaal zou zijn om op te nemen in de onderstaande methoden.

Een genenboom G op verzameling X van soorten is een boom met wortels (niet noodzakelijk binair) waarvan de bladeren zijn gelabeld (niet noodzakelijk bijectief) met X . Om de bladeren te onderscheiden die zijn gelabeld met hetzelfde element van X , voegen we subscripts toe aan de bladlabels. Figuur 1 toont een genenboom G begin x = <x, Y, Z> van soorten, waarbij vier allelen van soorten worden bemonsterd x (gelabeld x1, . x4), worden drie allelen bemonsterd uit soorten Y (gelabeld ja1, . ja3), en twee allelen worden bemonsterd uit soorten Z (gelabeld z1 en z2). In dit werk laten we in het bijzonder toe dat een genenboom nul allelen heeft die van sommige soorten zijn bemonsterd.

Genbomen en gewortelde triples. Een genenboom G op drie soorten X, Y , en Z, waarbij per soort meerdere allelen worden bemonsterd. Twee geïnduceerde triples door de genenboom en hun toewijzing aan de soortnamen worden getoond.

Een geroote triple (van nu af aan schrijven we gewoon "triple", aangezien we alleen te maken hebben met geroote topologieën) is een geroote boom met drie bladeren. Als de triple binair is, schrijven we xy|z om aan te geven dat de drievoudige put x en ja dichter bij elkaar dan een van hen z. Als de triple niet-binair is, dan is het xyz. We duiden aan met g|<x,y,z> de triple in de genenboom G geïnduceerd door de bladaanzet te beperken tot de drie gelabelde bladeren x, ja, en z. Figuur 1 toont de twee triples geïnduceerd door <x1, ja1, z1> en <x1, ja1, z2>. Ten slotte, om de bladlabels in de genenboom te koppelen aan hun overeenkomstige taxa in het fylogenetische netwerk, introduceren we functie φ die een allellabel in de genenboom toewijst aan het overeenkomstige taxon in het netwerk. In figuur 1 bijvoorbeeld φ(z1) = φ(z2) = Z. Verder gebruiken we φ(g|<x,y,z>) om de geïnduceerde tripel aan te duiden met zijn bladlabels vervangen door de taxa-namen (van soorten). Figuur 1 illustreert: φ.

Pseudo-waarschijnlijkheid van een soortennetwerk

Laat X een verzameling taxa (soorten) zijn, t = XY |Z wees een binair verdrievoudigen met X, Y, Z X , en G een genenboom zijn op X . We duiden aan met een(g, X), voor xx, de verzameling allelen van x dat label vertrekt van G. In figuur 1 bijvoorbeeld een(g, X) = <x1, x2, x3, x4>. wij definiëren ρ(t, g) om het aantal keren te zijn t wordt veroorzaakt door G (wanneer de bladlabels worden toegewezen aan X met behulp van de functie φ) genormaliseerd door het aantal keren dat een triple aan is gezet X, Y , en Z wordt veroorzaakt door G. Het is duidelijk dat als er maximaal één allel per soort wordt bemonsterd in G, dan wordt elke triple ofwel niet geïnduceerd door de genenboom of één keer geïnduceerd. Aangezien we echter meerdere allelen per soort toestaan, is dit mogelijk niet het geval. Merk op dat terwijl t binair is, kan het zo zijn dat g|<x l ,y J ,z k> is niet-binair. Aangezien er drie manieren zijn om een ​​niet-binair drietal op te lossen, is een niet-binair drietal g|<x l ,y J ,z k> draagt ​​1 . bij/3 tot de waarde van ρ(t, g). Rekening houdend met deze twee problemen, ρ(t, g) voor t = XY |Z gelijk aan

waarbij I de indicatorfunctie is gedefinieerd door I e = 1 wanneer e is waar en ik e = 0 wanneer e is fout. Voor een setje G van genenbomen definiëren we ρ ( t , G ) = ∑ g ∈ G ρ ( t , g ) . Als de noemer in Vgl. (1) is gelijk aan nul, we stellen ρ(t, g) = 0.

Gegeven een set G van genenbomen, de drie binaire triples t1 = XY|Z, t2 = XZ|Y, en t3 = YZ|X op een verzameling < X , Y , Z >⊆ X hebben een multinomiale verdeling gegeven door

waar P (t|Ψ,) is de kans op geroot triple t gegeven netwerk Ψ en overervingskansen γ [22, 21].

Ten slotte, de pseudo-waarschijnlijkheid van het fylogenetische netwerk Ψ en overervingskansen γ een set gegeven G van genenbomen wordt gegeven door

A maximum pseudo-likelihood approach seeks Ψ en γ 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 Ψ en γ .

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 , laten R(t) be the set of triples induced by tree t . Then no tree T′ exists such that T ≠ T′ en 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 en3 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, en B|CD. Branch lengths and inheritance probabilities are shown in blue and red, respectively, for Ψ1 en2.

A concrete example of these settings is:

network Ψ1: B1 = 1, B2 = 1, B3 = 2, B4 = 1, B5 = 0, α = 0.1

network Ψ2: ik1 = 1.841435, ik2 = 1.951019, ik3 = 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 Ψ ∗ en γ ∗

Given a set of gene trees G, Ψ en γ 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. (EEN) 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 NSxy. 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 SI-tekst.

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. Bijvoorbeeld, 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. (EEN) Progression over time from population splitting, divergence in isolation, and secondary contact, with alternate outcomes: (l) tension zone, (ii) merging, and (iii) overriding of expanding population (blue) over the resident population with introgression from yellow→blue for some genes. (B en 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 FNS (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? - 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.


Methoden:

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

Laten 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, namelijk 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,ja) ∈ E(N) | <x,ja> ⊆ V(N) ∖ V(S)> for a subnetwork S van N.

Voor jij,vV(N), jij is een parent van v en v is een kind van jij if (jij,v) ∈ E(N). We gebruiken C(R) to denote the unique child of RR(N). If there is a direct path from jij tot v, v heet a descendant van jij.

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 onderstaand R, we use N−[ R] N+ to denote the subnetwork obtained by replacing [ R] N met zodat 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 binair 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 verminderd.

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, een TROS 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 over 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 over x and a phylogenetic tree t over x, we say that t is displayed in N indien 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 met wie which has two incident edges (jij,met wie) en (met wie,v), we merge the two edges into one edge (jij,v).

A cluster Bx is een soft cluster in N if there is a tree t displayed in N zoals dat 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. We gebruiken S C(N) to denote the set of soft clusters in N.

Gegeven Bx and a network t Aan x, the CCP asks whether B is a soft cluster in N [10], which is formulated as below:

Voorbeeld: A phylogenetic network N over a set of taxa x en Bx.

Vraag: Is BS C(N)?

Laten N 1 en 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 NS 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) en 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 zichtbaar on a leaf if it lies on all the paths from ρ(N) tot . If a node RR(N) ∪ <ρ(N)> is visible on a leaf , its tree component C R is zichtbaar Aan ook. 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 blootgesteld if it contains only one leaf or if all the tree components right below it are trivial.

Blijkbaar, 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, en C r6, waar C r6=. C r7 and r7 are right below C r5, C r2, en C R. C R is visible on all the leaves. C r1 en 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 en N 2 are derived from N, which contain fewer nodes than N. The algorithm is then called on both instances ( C,N 1 ) en ( 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.

Laten N be a network over x en Bx, respectievelijk. We examine a non-trivial exposed tree component C R van 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. We gebruiken l R(C R) en C R(C R) to denote the sets of inner- and cross- reticulation nodes, respectively. For example, in Fig. 1, l R(C r5)= en C R(C r5)=.

We gebruiken 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 , ekst >) when assuming B=.

Stel dat 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 en B.

Let (ar =X ackslash B) . Indien L RB and ( L_ cap ar eq emptyset ) , B must be a soft cluster of a node in C R indien 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:

Sinds L RB and (check _ subseteq B) , (hat subseteq B) . If (hat = B) , B is a soft cluster of ρ(C R) in N een. Otherwise, if (hat subset B) , we set:

Indien 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 en 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

Aannemen dat C R is a non-trivial, exposed and visible tree component in a network N over the taxa set x, en dat Bx. Laten L R, (hat ) , B ′ , N een′, and N B′ be defined above.

(i) If (hat subset B) , B is a soft cluster in N als en alleen als B ′ is a soft cluster in N een′.

(ii) Indien B is not a soft cluster of a node in C R en L RB= , B is a soft cluster in N als en alleen als B is a soft cluster in N B′.

Stel dat C R is not visible. Indien 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. Indien C R=<C(R)> and (c(r)=r'phantom !>) , then at least one parent of R ′ is not R. We can now define: