Det finns många gemensamma polymorfismer mellan A. thaliana och C. rubella
I en population av 80 A. thaliana accessioner fanns det 4 902 039 SNP:er (av 119 146 348 platser), bland vilka 2 044 731 hade en mindre allelfrekvens (MAF) av > 0,05. I C. rubella-populationen identifierade vi 2 149 643 SNP:er (av 134 834 574 platser) genom att kalla SNP:er från 22 C. rubella-accessioner (Additional file 1: Table S1, inklusive 21 publicerade accessioner och en accession som sekvenserats i den här studien) mot C. rubella-referensgenomet, av vilka 1 240 547 hade en MAF > 0,05. För att identifiera gemensamma polymorfismer mellan de två arterna, definierade som samma allelpar på en viss ortolog plats, konstruerade vi först uppsättningen ortologa genpar mellan de två arterna. För att garantera att de ortologa generna är bevarade inkluderade vi förutom referensgenomerna för A. thaliana och C. rubella även Arabidopsis lyrata , en släkting till A. thaliana. Vi fick 16 047 ortologa genpar och tog bort 33 som hade tandemduplikationer i någon av de tre referenserna och fick slutligen totalt 16 014 ortologa genpar mellan A. thaliana och C. rubella för vidare analys.
Den geniska regionen för de 16 014 ortologa generna i A. thaliana sträckte sig över 39 275 210 bp och på samma sätt sträckte den sig över 40 936 262 bp i C. rubella. Dessa regioner innehöll 3 889 495 fasta skillnader och denna höga andel (~ 10 %) är förenlig med den långa divergenstiden (~ 8 MYA) mellan de två arterna. I dessa regioner hittade vi 1 122 845 bi-alleliska platser (426 123 med MAF > 0,05) i A. thaliana och 452 116 bi-alleliska platser (279 780 med MAF > 0,05) i C. rubella. Bland dessa polymorfa platser var 19 732 ortologiska platser polymorfa i båda arterna, varav 8535 delade samma allelpar (shared SNP ) (Additional file 1: Table S2).
Vid jämförelse med sekvenser i icke-kodande regioner är sekvenser i kodande regioner mer bevarade och ger robusta anpassningar mellan de två starkt divergerande arterna; därför fokuserade vi först på shSNPs i kodande regioner. MAF > 0,05 krävdes i båda arterna för att garantera SNP-tillförlitlighet och ta hänsyn till det förväntade överskottet av alleler med intermediära frekvenser för platser som är föremål för långsiktigt balanserande urval. Vi hittade 1503 shSNP i de kodande regionerna av 1007 gener.
För att undvika genotypnings- och kartläggningsfel tillämpades ytterligare filtrering på de 1503 shSNP. Filtreringen tillämpades endast på C. rubella SNP-data, eftersom vi hämtade SNP-matrisen för A. thaliana. För att undvika falska SNP:er till följd av duplikationer i genomet bedömde vi kartläggningsbarheten för varje 50-bp-region i C. rubella och behöll endast platser som låg i unikt kartläggningsbara regioner för efterföljande analys. Detta lämnade endast 580 platser kvar. Slutligen, efter att ha tagit bort platser av låg kvalitet som markerats av SNP-kallningsverktyget, fick vi 546 tillförlitliga gemensamma kodande SNP:er i 433 gener. Detaljer om filtreringsprocessen finns i avsnittet ”Metoder” och en bild av processen visas i fig. 2.
Pipeline för SNP-filtreringsprocessen för att identifiera TSP-kandidatplatser
De två arternas demografiska historia
Detektering av riktiga TSP-signaler från de rikligt förekommande delade polymorfismerna är beroende av en fullständig förståelse av de två arternas demografiska historia. Det gemensamma frekvensspektrumet (Joint Site Frequency Spectrum, SFS) har använts i stor utsträckning för att studera den demografiska historien hos olika organismer . Därför extraherade vi först de fyrfaldigt degenererade platserna från anpassningarna av referensgenomerna för A. thaliana och C. rubella på de 16 014 ortologerna. Slutligen fick vi fram 2 011 573 platser för den demografiska analysen (se ”Metoder” för mer information).
Simuleringar av koalescens utfördes sedan med hjälp av fastsimcoal2 enligt en grundmodell utan genflöde (M1, fig. 3) och en modell som innehåller gammalt genflöde mellan de två släktena (M2, fig. 3). Vi tog endast hänsyn till gammalt genflöde mellan de två arterna, eftersom det är högst osannolikt att arter som tillhör olika släkten och som har olika antal kromosomer (fem respektive åtta) har introgression nyligen. Dessutom är A. thaliana i båda släktena den enda arten med fem i stället för åtta kromosomer ; vi begränsade därför det gamla genflödet innan A. thaliana separerades från resten av Arabidopsis-släktet. I varje modell fastställde vi divergenstiden för de två släktena till 8 MYA , vilket motsvarar 8 miljoner generationer sedan, och antog en spontan mutationshastighet på 7 × 10-9 per bp per generation . Vi har övervägt olika populationsstorlekar för båda arterna baserat på övergångshändelserna från deras respektive stamfäder; A. thaliana genomgick en populationsminskning efter att den divergerade från resten av Arabidopsis-släktet omkring 6 MYA och C. rubella upplevde en mycket nyligen inträffad flaskhals i samband med artbildningen från C. grandiflora . Vi använde koalescenssimuleringar med tillämpning av den sammansatta sannolikhetsmetoden som implementerats i fastsimcoal2 för att anpassa båda modellerna till den gemensamma SFS för de två arterna som beräknats från de extraherade 2 011 573 fyrafaldigt degenererade platserna över artgränserna. Vi jämförde de två modellerna med hjälp av Akaikes informationskriterium (AIC) och Akaikes bevisvikt (w), som i Excoffier et al. . Modellen utan gammalt genflöde (M1) passade något bättre (Max EstLhood: -682010 vs -682028), med ett lägre AIC och högre vikt än för den andra modellen (fig. 3, Additional file 2: Table S3). Dessutom tyder de två nära sannolikheterna på att effekten av det förfäderliga genflödet borde ha utplånats över den långa tidsskalan och bidrar föga till modellens kvalitet.
Skattningar av demografiska parametrar för två modeller av de två arternas divergens
Under modell M1 är den nuvarande N e av A. thaliana var ~ 519 000 med ett 95 % konfidensintervall (CI) på 486 368-527 574, från en stor förfäderspopulation (~ 2 230 000, 95 % CI = 1 085 330-4 876 051) innan den separerades från resten av Arabidopsis-släktet vid ~ 5,84 MYA (95 % CI = 5,27-6,70). C. rubella utvecklades ~ 0,40 MYA (95 % CI = 321 998-500 317) från en förfäderspopulation med en stor N e på ~ 4 037 000 (95 % CI = 2 076 868-5 165 614) och en nuvarande N e på ~ 129 000 (95 % CI = 126 383-157 779). De två släktena har divergerat från en stampopulation med N e = ~ 4 930 000 (95 % KI = 4 560 931-4 969 696). Enligt modell M2 med genflöde erhölls liknande parameteruppskattningar, med undantag för ett större ursprungligt N e för Arabidopsis-släktet (~ 3 270 000, 95 % KI = 797 016-4 342 346) och ett mindre N e för Capsella-släktet (~ 1 972 000, 95 % KI = 2 126 346-6 248 003). Ett starkare genflöde uppskattades från Capsella till Arabidopsis än i motsatt riktning (migrationshastighet per generation; 1 × 10-8, 95 % CI = 4,0 × 10-15-1,1 × 10-6 jämfört med 7 × 10-14, 95 % CI = 5,7 × 10-15-6,1 × 10-5), även om båda var svaga (se tilläggsfil 2: tabell S3 för detaljer).
Trans-species polymorfismer mellan de två arterna måste vara under balanserande urval
Trans-species polymorfismer kan vara neutrala och dess sannolikhet kan approximeras givet specifika demografiska parametrar. I likhet med en studie av TSPs hos människor och schimpanser , under neutral evolution, var delade polymorfismer identiska genom härstamning i vårt system endast om: (1) minst två A. thaliana-linjer och två C. rubella-linjer inte samlades före uppdelningen mellan A. thaliana och C. rubella, och (2) linjer som bär på samma allel samlades före linjer som bär på olika alleler. Denna sannolikhet bestäms huvudsakligen av villkor (1) och kan approximeras genom följande baserat på koalescensteorin :
där T avser divergenstiden för de två släktena och N A/N C avser populationsstorlekarna hos A. thaliana/C. rubella. Enligt våra uppskattningar enligt modell M1, med beaktande av förändringar i populationsstorlek, är denna sannolikhet för identitet genom härstamning i storleksordningen 10-9. Med tanke på att vi har < 39 275 210 anpassade platser mellan de två arterna i den geniska regionen, förväntar vi oss att det totala antalet neutrala TSP:er blir < 1 enbart genom genetisk drift.
Vi antog slumpmässig parning i vår modell; båda arterna är dock självbildande och det finns troligen en populationsstruktur inom arterna. Trots detta bör nyligen inträffade demografiska händelser ha relativt liten effekt, eftersom vi kräver djupa koalescenshändelser av slump hos båda arterna i samma region av genomet . Som illustrerats i tidigare studier bör även djup populationsstruktur inom moderna människor ha minimal effekt på sannolikheten. I den här studien har båda arterna en historia av övervägande korsning. A. thaliana övergick från utkorsning till självkorsning för bara en miljon år sedan, medan C. rubella övergick från utkorsning till självkorsning för mycket senare tid sedan. Även när det gäller självsvältande arter är utkorsningsgraden i lokala populationer så hög som 14,5 % . Därför är det osannolikt att populationsstrukturer, om de finns, kommer att bestå över en lång tidsskala och dess inverkan på sannolikheten kan därför ignoreras.
Identifiering av artöverskridande polymorfismer under balanserande urval
TSP:er kan särskiljas från neutrala mutationer eftersom regioner som är föremål för långvarigt balanserande urval grupperar sig per allel, snarare än per art . Därför fokuserade vi nästa gång på de 433 kandidatgenerna med tillförlitliga delade SNP:er i den kodande regionen och undersökte haplotyperna som täcker varje delad bi-allelisk SNP med MAF > 0,05 i de geniska regionerna.
För att uppskatta längden på varje segment som bär på en signal av TSP:er använde vi en formel som tidigare härletts och som till stor del är beroende av rekombinationshastigheten. Ur koalescenssynpunkt bryts ett sådant segment inte upp av rekombination förrän alla linjer från samma allelklass koalescerar till sin senaste gemensamma förfader i förfäderspopulationen . Om man antar en rekombinationshastighet på 3,6 cM/Mb för båda arterna var segmentets längd extremt kort, dvs. teoretiskt sett bara några baspar. Med tanke på att båda arterna nyligen uppstod från sina respektive korsande stamfäder och att den effektiva rekombinationshastigheten kan ha varit mycket högre tidigare, kan den förväntade längden vara ännu kortare. Denna uppskattning antyder, under de neutrala omständigheterna i vårt system, att det är mycket svårt att upptäcka något segment utan ett avbrott i rekombinationen. När det finns ett balanserande urval kan dock urvalet undertrycka rekombination i det omgivande området . Därför bör segmentlängden vara längre än den teoretiska uppskattningen enligt en neutral modell. Vi skannade därmed den geniska regionen med en fönsterstorlek på 100 bp och en stegstorlek på 1 bp.
I de 433 kandidatgenerna upptäckte vi 975 delade bi-alleliska SNPs (inklusive både exoniska och introniska SNPs med MAF > 0,05). I likhet med tidigare studier letade vi sedan efter fönster som täcker minst två av de 975 SNP:erna som befinner sig i stark kopplingsobalans (r 2 > 0,5) hos båda arterna bland de kvalificerade fönstren (anpassade till minst 95 % av längden; se ”Metoder” för detaljer) för att identifiera alleliska träd. Dessa begränsningar kan kraftigt minska antalet falska positiva resultat och ge allelträd, om de finns, med hög upplösning. Slutligen identifierade vi fönster från fem gener, AT1G35220, AT2G16570, AT4G29360, AT5G38460 och AT5G44000, som omfattar tio platser, som kandidat-TSP:er under långvarigt balanseringsurval (Additional file 3: Figur S1). Ingen av de fem ortologa generna som vi hittade här är korrelerade med variation i kopianummer (CNV) och alla har endast en träff när vi jämförde dem mot referenserna för de två arterna, respektive (se ”Metoder” för detaljer).
För att verifiera de identifierade regionerna bestämde vi först alla haplotyper i de identifierade regionerna från varje population och resequencerade representativa accessioner för varje haplotyp (se Additional file 1: Tabell S4 för primers). Som förväntat validerades alla TSP-kandidatplatser i de fem generna och sekvenserna hos de två arterna i kandidatregionerna klustrades efter allel, snarare än efter art (fig. 4). I genen AT1G35220 befann sig de två TSP-kandidatplatserna i fullständig kopplingsdifferens i en intronisk region; denna region kan vara måltavla för balanserande urval eller kopplad till en oupptäckt kodande TSP-plats.
Alla kandidatregioner i de fem generna ger ett allelträd, snarare än ett artträd
Och även om haplotyperna för varje region klustrade efter allel, snarare än art, upptäcktes sällan delning av haplotyper mellan de två arterna, utom i AT2G16570 (Col-0 delade sin haplotyp med flera C. rubella-tillgångar; fig. 4). Detta är inte förvånande med tanke på den långa divergenstiden; omfattande haplotypdelning uppträder vanligtvis på en mycket kortare tidsskala och induceras av händelser som nyligen inträffad introgression mellan närbesläktade arter.
Neutrala simuleringsstudier validerar de fem kandidatgenerna
För att se om de observerade fönstren skulle kunna genereras slumpmässigt under neutral evolution, vilket skulle resultera i falska positiva resultat, körde vi ytterligare simuleringar baserade på de uppskattade demografiska parametrarna med hjälp av fastsimcoal2 (Additional file 4: Text S1). Förutom neutrala återkommande mutationer kan genflöde också resultera i delade SNPs. Följaktligen körde vi simuleringar enligt både modell M1 (utan genflöde) och M2 (med gammalt genflöde), även om vår demografiska analys visade att M1 stämde något bättre med data. I båda simuleringarna tog vi hänsyn till heterogenitet i mutationsfrekvenser för olika klasser av mutationer, särskilt den högre mutationsfrekvensen på CpG-platser, vilket kan resultera i falskt positiva resultat (Additional file 1: Table S5, Additional file 4: Text S1). Med hjälp av fastsimcoal2 genererade vi 1 000 000 neutrala segment på 100 bp under varje modell och letade efter dem med två eller fler delade SNP:er och kluster per allel som vi sökte efter TSP:er.
För båda modellerna gav ingen av de 1 000 000 körningarna upphov till ett fönster som uppfyllde våra kriterier (Additional file 1: Table S6). Trots förekomsten av neutrala delade SNP:er gav inget simulerat fönster upphov till ett allelträd, eftersom alla fönster med delade SNP:er åtföljdes av mycket mer fasta skillnader mellan de två arterna, vilket antyder högre divergensnivåer än mångfald. Detta resultat tyder på att dessa simulerade neutrala delade SNP:er är återkommande mutationer snarare än TSP:er, och ännu viktigare är att de fem gener vi hittade inte är förenliga med neutral evolution och därmed visade sig vara riktiga TSP:er under balanserande urval. De slutliga TSP-platserna och generna anges i tabell 1. Tillsammans med den tidigare nämnda demografiska studien innebär våra resultat dessutom att även om forntida genflöde förekom, skulle TSP:erna under neutral evolution gå förlorade genom drift i detta system.
Egenskaper hos generna under balanserande urval
Vi beräknade därefter nukleotiddiversiteten (π) för alla TSP-regioner i de fem generna hos varje art och använde de simulerade neutrala sekvenserna under M1 för att bestämma bakgrundsdiversitetsnivåer. Alla regioner i de fem generna uppvisade signifikant högre π-värden än bakgrundsnivåerna i både C. rubella och A. thaliana (Wilcoxon-Mann-Whitney-test, FDR-korrigerat P < 0,05, Tabell 2, Additional file 3: Figur S2A), utom AT5G38460 i A. thaliana. Dessutom visade allelerna i dessa gener en tendens till intermediära frekvenser (Wilcoxon-Mann-Whitney-test, P = 0,0752/0,03474 för A. thaliana/C. rubella; Additional file 3: Figur S2B). En intermediär frekvens är dock en indikation på balanserande urval, men inte ett definitivt bevis, eftersom allelfrekvensfördelningen för platser som är kopplade till en balanserad polymorfism förväntas uppvisa en förskjutning mot frekvensjämvikten, som kan vara vid vilken allelfrekvens som helst.
En av de fem generna under långvarigt balanserande urval i den här studien, AT1G35220, har en okänd funktion, men uppvisar proteinfosforylering under etylenbehandling . AT2G16570 är bland annat ett nyckelenzym i purinnukleotidbiosyntesen och är viktig för celldelning, kloroplastbiogenes och fröspridning. AT4G29360 är ett protein från O-glykosylhydrolasets 17:e familj, som är involverat i försvarsreaktioner; AT5G38460 är ett glykosyltransferas som katalyserar överföringen av en glykosylgrupp från en förening (donator) till en annan (acceptor) och är involverad i olika funktioner, inklusive biotisk stress ; AT5G44000 är ett glutation S-transferas, som vanligtvis är involverat i svaret på abiotisk och biotisk stress . Uppenbarligen är dessa gener potentiellt involverade i responsen på biotisk eller abiotisk stress (AT4G29360, AT5G38460 och AT5G44000) eller grundläggande biokemiska funktioner (AT2G16570).
Som väntat var de gener som var föremål för balanserande urval funktionellt viktiga och alla homologer till de fem generna fanns redan i de gröna växternas senaste gemensamma förfader. Som framgår av tabell S7 (Additional file 1: Table S7) kan homologer (antingen ortologer eller paraloger) hittas till och med i den mest basala arten av gröna växter, Chlamydomonas reinhardtii, för alla de fem generna, utom AT4G29360, som kan spåras tillbaka till Physcomitrella patens.
Däremot utmärkte sig inte loci som allmänt accepteras vara under balanserande urval, såsom S-locus eller R-gener , i denna studie. Detta är förväntat, eftersom dessa loci är alltför variabla för att kunna identifieras baserat på korta läsningar. Till exempel är R-gener för dynamiska för att kallas SNPs; S-locus existerar inte i den senaste annoteringen av Arabidopsisgenomet och endast en S-locus haplotyp upprätthålls i C. rubella sedan övergången från outcrossing till selfing och upplösningen av självinkompatibiliteten . Dessutom är S-locus inte längre föremål för balanserande urval, eftersom båda arterna nu är självutnyttjande. Däremot har de gener vi identifierat här, även om de är gamla, inte studerats utförligt och kan ge en inblick i vilka typer av gener som står under balanserande urval.
Balanserande urval bidrog till anpassning till divergerande livsmiljöer
För att se om de allelvarianter som står under långvarigt balanserande urval är förknippade med ekologisk diversifiering undersökte vi divergensen med avseende på 48 ekologiska faktorer (Additional file 5: Table S8A). På grund av brist på GPS-information och den lilla provstorleken för C. rubella var denna analys endast möjlig för A. thaliana-proverna. Populationsstrukturen är vanligtvis starkt korrelerad med ekologisk diversifiering och kan därför förvränga våra resultat. Vi kontrollerade först om någon TSP-plats var korrelerad med populationsstrukturen i A. thaliana-proverna, även om en sådan struktur inte påverkar sannolikheten att observera artträdet för A. thaliana och C. rubella. Med hjälp av ADMIXTURE fann vi att de 80 A. thaliana-proverna kan klassificeras i två grupper (Additional file 3: Figur S3; Additional file 6: Tabell S9) och att endast de alleliska klassificeringarna av de två platserna från genen AT5G38460 är signifikant korrelerade med populationsstrukturen (chi-två-test, FDR-korrigerad P < 0,05,; Additional file 1: Tabell S10). Vi uteslöt därmed AT5G38460 från efterföljande ekologiska analyser.
För att få en grundlig förståelse för ekologisk divergens använde vi 1135 nyligen publicerade genomer av A. thaliana . Först tillämpade vi en ”gallringsprocess” för att garantera att varje prov var mycket representativt för sin naturliga livsmiljö, vilket lämnade 584 prov (se ”Metoder”). För det andra klassificerade vi för varje gen de 584 accessionerna av A. thaliana i två grupper baserat på de fasade haplotyperna för de två TSP-platserna (Additional file 5: Table S8B, C, vissa prover togs bort eftersom de inte kunde fasas). Vi utvärderade sedan divergensen mellan de två grupperna av accessioner med avseende på de 48 ekologiska faktorerna för var och en av de fyra generna. Intressant nog var alla dessa fyra gener förknippade med divergens för vissa specifika ekologiska parametrar. Särskilt AT1G35220 och AT4G29360 uppvisade signifikant divergens med avseende på de flesta av de temperaturrelaterade ekologiska faktorerna (Additional file 5: Table S8 A, Wilcoxon-Mann-Whitney-test, FDR-korrigerat P < 0,05).
Vi modellerade därefter de ekologiska nischerna för alla fyra generna. Uppenbarligen uppvisade de två grupperna av prover för varje gen, vilket indikeras av Warrens I-statistik som mäter nischlikhet , betydligt lägre observerad nischidentitet än 100 slumpmässiga permutationer (t-test med ett urval, FDR-korrigerat P < 0,01; Fig. 5a, Additional file 5: Table S8 D). Med andra ord uppvisar de två allelgrupperna av prover betydande nischdivergens. Dessutom var proverna av varje allelisk typ för varje gen utspridda, i stället för att vara isolerade i ett litet lokalt område (Additional file 3: Figur S4). Dessa resultat tyder på att alla dessa loci är korrelerade med anpassning.
Ekologisk och uttrycksdivergens. a Betydande ekologisk divergens mellan de två typerna av prover för var och en av de fyra generna, vilket framgår av den observerade I-poängen (I O) och simulerade I-poängen (I S). b Expressionsdivergens för genen AT5G44000. c Till vänster: Modellering av nischen med hög sannolikhet (≥ 0,5) för de två typerna av prover för AT5G44000. Höger: Signifikansresultat under olika permutationsstrategier (för nischer med sannolikhet ≥ 0,5; I O = 0,673, 100 permutationer)
Vi undersökte också uttrycksdifferentiering för de fyra generna mellan de två motsvarande grupperna baserat på de fasade haplotyperna vid de två TSP-platserna genom att välja 84 publicerade bladvävsextraherade transkriptomer från A. thaliana (ett prov sekvenserades för varje accession och uttrycksnivån mättes som fragment per kilobas av exon per miljon kartlagda fragment ) som vår tidigare studie . En gen, AT5G44000, uppvisade en signifikant uttrycksskillnad (Wilcoxon-Mann-Whitney-test, FDR-korrigerat P < 0,05, fig. 5b) mellan de två haplotypgrupperna.
Vi utförde därför en djupgående nischmodellering av AT5G44000 (fig. 5c) och undersökte diversifieringen av de två grupperna av prover (503 vs 75). Vi jämförde först nischidentiteten mellan de två haplotypgrupperna av AT5G44000 genom att begränsa vår analys till nischer med hög sannolikhet (≥ 0,5) och fick liknande resultat (Fig. 5c, Additional file 5: Table S8 D). För att se om den obalanserade provstorleken kunde påverka resultaten använde vi en annan permutationsstrategi genom att begränsa analysen till samma provstorlek (75) för båda uppsättningarna i varje upprepning (med sannolikhet > 0,5). Som visas i figur 5c visade det observerade I-värdet (0,673) ingen signifikant skillnad när permutationen utfördes för de verkliga provgrupperna (simulering 1) (t-test med ett urval, P = 0,166), vilket tyder på att det observerade värdet var tillförlitligt, oberoende av skillnaden i provstorlek. När de två verkliga grupperna blandades och två slumpmässiga grupper av verklig storlek valdes ut (simulering 2) eller två slumpmässiga grupper av samma storlek (75) valdes ut (simulering 3), var skillnaden mellan det observerade värdet och permutationerna återigen signifikant (t-test med ett urval, P = 1,9 × 10-75 för simulering 2 och P = 2,6 × 10-75 för simulering 3). Dessa resultat innebär att de två funktionellt differentierade haplotypgrupperna av AT5G44000 har anpassat sig till olika ekologiska livsmiljöer.