Sampling

Kotieläiminä pidetyistä lampaista (O. aries) ja vuohista (C. hircus) otettiin näytteet Iranissa (IROA- ja IRCH-ryhmät) ja Marokossa (MOOA- ja MOCH-ryhmät), ja näytteitä kerättiin kaikkiaan 20 eläimestä ryhmää kohden (Supplementary Fig. 6). Näytteet kerättiin tammikuun 2008 ja maaliskuun 2012 välisenä aikana Marokon pohjoisosassa ja elokuun 2011 ja heinäkuun 2012 välisenä aikana Luoteis-Iranissa eurooppalaisen Nextgen-hankkeen puitteissa (avustussopimus nro 244356) Euroopan unionin direktiivin 86/609/ETY eettisten sääntöjen mukaisesti. Korvaleikkeet kerättiin satunnaisesti valittujen eläinten korvan distaalisesta osasta, ja niitä säilytettiin välittömästi 96-prosenttisessa etanolissa yhden päivän ajan ennen kuin ne siirrettiin silikageelihelmiin DNA:n uuttamista varten.

Luonnonvaraisista lajeista aasialaisesta muflonista (O. orientalis) ja bezoar-mäyräkoirasta (C. aegagrus) otettiin näytteet Luoteis-Iranissa kotieläinjalostuksen kehtoon kuuluvista lajeista21,22. Kolmetoista aasialaista muflonia ja 18 bezoar-mäyräkoiran kudosta (IROO- ja IRCA-ryhmät, täydentävä kuva 6) kerättiin joko vankeudessa pidetyistä tai hiljattain metsästetyistä eläimistä sekä Iranin ympäristöministeriön saatavilla olevista pakastetuista näytteistä. Tämän yksilöpohjaisen näytteenottomenetelmän tarkoituksena on minimoida mahdolliset vääristymät välttämällä paikallisten vaikutusten yliedustus (esim,

Lisätiedot

Lisäksi koottiin maailmanlaajuinen rotupaneeli lampaille ja vuohille (wpOA ja wpCH). wpOA sisälsi 20 WGS-näytettä 12-kertaisella kattavuudella, jotka edustivat 20:tä erilaista maailmanlaajuista rotua. wpOA koostui 14:stä WGS-näytteestä, jotka oli sekvensoitu 12-kertaisella kattavuudella ja jotka edustivat 9:ää eurooppalaista yksilöä, ts, 2 ranskalaista alppilammasnäytettä ja 2 ranskalaista saaninäytettä, jotka INRA on sekvensoinut, 5 italialaista saaninäytettä, jotka Parco Tecnologico Padano on toimittanut, ja 5 australialaista yksilöä, ts, 2 Boer-, 2 Rangeland- ja 1 Cashmere-näytettä, jotka CSIRO toimitti (Supplementary Data 5).

WGS-datan tuottaminen

Genominen DNA saatiin onnistuneesti uutettua kaikista kudosnäytteistä Macherey Nagel NucleoSpin 96 Tissue kitillä valmistajan protokollaa soveltaen. Kudosnäytteenotto suoritettiin MN-neliökuoppalohkoihin, jotta saatiin 25 mg fragmentteja näytettä kohti. Valmistettiin kolme ja puoli MN-neliö-96-lohkoa, ja uutto suoritettiin Tecan Freedom Evo Liquid handler -laitteella valmistajan protokollan mukaisesti. Esilysointivaiheessa näytteet homogenoitiin 180 µl:lla T1-puskuria ja 25 µl:lla proteinaasi K:ta yön yli 56 °C:ssa. Sitoutumisolosuhteiden säätämiseksi lisättiin 200 µl BQ1-puskuria ja näytelevyä inkuboitiin 1 h 70 °C:ssa; tämän jälkeen lisättiin 200 µl 100-prosenttista etanolia. Lysaatit siirrettiin Nucleospin Tissue -sidontalevylle ja läpivirtauksen poistamiseksi käytettiin tyhjiötä (-0,2 bar, 5 min). Kolme pesuvaihetta tehtiin vastaavasti BW- ja B5-puskureilla, ja läpivirtaus poistettiin jälleen tyhjiöllä. Ennen genomisen DNA:n eluointia Nucleospin Tissue binding plate -silikamembraani kuivattiin tyhjiössä vähintään -0,6 baarissa 10 minuutin ajan. Eluointivaihe suoritettiin 100 µl:lla esilämmitettyä BE-puskuria (70 °C) ja sentrifugointivaiheella 3700 rcf:n kierrosluvulla 5 minuutin ajan 96-PCR-kuopissa. Genomista DNA:ta säilytettiin 4 °C:ssa jäätymissulatuksen välttämiseksi ja sen konsentraatio (ng/μl:nä) testattiin Picogreen-menetelmällä ja Nanodrop-laitteella.

Kokonaiset genomit sekvensoitiin uudelleen 500 ng:sta genomista DNA:ta, jotka leikattiin 150-700 bp:n alueelle Covaris® E210 -laitteistolla kutakin näytettä kohden ja käytettiin Illumina®-kirjastojen valmistukseen puoliautomatisoidun protokollan avulla. Loppujen korjaus, A-tailing ja Illumina®-yhteensopivien adapterien (BioScientific) ligointi suoritettiin SPRIWorks-kirjastonvalmistusjärjestelmällä ja SPRI TE -laitteella (Beckmann Coulter) valmistajan protokollan mukaisesti. Useimpien fragmenttien talteenottamiseksi käytettiin 300-600 bp:n kokovalintaa. DNA-fragmentit monistettiin 12 PCR-syklin avulla käyttäen Platinum Pfx Taq Polymerase Kit -pakettia (Life® Technologies) ja Illumina®-adapterispesifisiä alukkeita. Kirjastot puhdistettiin 0,8x AMPure XP -helmillä (Beckmann Coulter) ja analysoitiin Agilent 2100 -bioanalysaattorilla (Agilent® Technologies) ja qPCR-kvantifioinnilla. Kirjastot sekvensoitiin käyttäen 100 emäspituuden lukukemiaa pareittaisessa virtaussolussa Illumina® HiSeq2000 -laitteella.

Illuminan pareittaiset lukemat kartoitettiin Ovis-lajin osalta lampaan referenssigenomiin (OAR v3.1, GenBank-kokoonpano GCA_000298735.146) ja Capra-lajin osalta vuohen referenssigenomiin (CHIR v1.0, GenBank-kokoonpano GCA_000317765.147) BWA-MEM48:n avulla. Kullekin yksilölle tuotettu BAM-tiedosto lajiteltiin Picard SortSam -ohjelmalla ja parannettiin käyttäen sekventiaalisesti Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator- ja GATK IndelRealigner49 -ohjelmia sekä SAMtools calmd50 -ohjelmia.

Varianttien löytämiseen käytettiin kolmea eri algoritmia: Samtools mpileup50, GATK UnifiedGenotyper51 ja Freebayes52. Varianttipaikat tunnistettiin itsenäisesti kullekin kuudesta ryhmästä käyttäen kutsualgoritmien moniotosmoodeja: (i) 162 näytettä MOOA:sta; (ii) 20 näytettä IROA:sta; (iii) 14 näytettä IROO:sta; (iv) 162 näytettä MOCH:sta; (v) 20 näytettä IRCH:sta; (vi) 19 näytettä IRCA:sta. Joidenkin ryhmien osalta useampien yksilöiden WGS oli saatavilla osana NextGen-hanketta (ks. edellä). Tässä tutkimuksessa käytetyt näytteet valittiin siten, että saatiin mahdollisuuksien mukaan tasapainoisia 20 yksilön ryhmiä. IRCA- ja IROO-ryhmien osalta lisänäytteitä saatiin myöhemmässä vaiheessa, ja ne lisättiin myöhempiä analyysejä varten. Lopullisen aineiston (Supplementary Data 5) saamiseksi poistettiin eläimet, joiden kohdistus- ja kutsumislaatu oli heikko.

Jokaisessa ryhmässä tehtiin kaksi peräkkäistä kierrosta varianttipaikkojen laadun suodatusta. Suodatusvaiheessa 1 yhdistettiin kolmesta algoritmista saadut kutsut ja suodatettiin samalla matalimman luotettavuuden omaavat kutsut pois. Varianttipaikka läpäisi suodatuksen, jos sitä kutsuttiin vähintään kahdella eri kutsualgoritmilla, joiden phred-variantin laatu oli >30. Vaihtoehtoinen alleeli paikassa hyväksyttiin, jos se kutsuttiin jollakin kutsuvalla algoritmilla ja genotyyppiluku oli >0. Suodatusvaiheessa 2 käytettiin GATK:n suorittamaa Variant Quality Score Recalibrationia. Ensin luotiin harjoitusjoukko korkeimman varmuuden omaavista varianttikohteista ryhmässä, jossa (i) kaikki kolme varianttikutsuohjelmaa kutsuvat kohdetta phred-asteikollisella varianttilaadulla >100, (ii) kohde on biallelinen, (iii) pienempien alleelien määrä on vähintään 3, kun lasketaan vain näytteet, joiden genotyyppi phred-asteikollisella laadulla >30. Koulutusjoukkoa käytettiin Gaussin mallin rakentamiseen GATK VariantRecalibrator -työkalulla käyttäen seuraavia UnifiedGenotyperin varianttiannotaatioita: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Gaussin mallia sovellettiin koko aineistoon, jolloin saatiin VQSLOD (log odds ratio of being a true variant). Kohteet suodatettiin, jos VQSLOD <cutoff-arvo. Cutoff-arvo asetettiin kullekin ryhmälle seuraavasti: VQSLOD-minimi = {harjoitusjoukon varianttien VQSLOD:n mediaaniarvo}-3 × {harjoitusjoukon varianttien VQSLOD:n absoluuttisen poikkeaman mediaani}. Siirtymä/transversio SNP-suhde viittasi siihen, että valittu rajauskriteeri antoi parhaan tasapainon selektiivisyyden ja herkkyyden välillä.

SNPs-kutsusarjat luotiin kuudelle Ovis- ja Capra-eläinryhmälle (eli iranilaisille ja marokkolaisille kotieläimille ja luonnonvaraisille eläimille kummallekin suvulle). Koska tässä tutkimuksessa tehdyt analyysit edellyttivät ryhmien välisiä vertailuja, luotiin genotyyppikutsusarjat yhdenmukaisella SNP-kohtien joukolla kaikille eläimille mistä tahansa ryhmästä. Yhdistimme kunkin suvun osalta sen kolmen ryhmän varianttikutsusivustot ja säilytimme vain bialleliset paikat, joista ei puuttunut tietoja. Genotyypit kutsuttiin uudelleen kussakin biallelisessa SNP-kohdassa kaikille kiinnostuksen kohteena oleville yksilöille GATK UnifiedGenotyper -ohjelmalla käyttäen vaihtoehtoa GENOTYPE_GIVEN_ALLELES. Tässä vaiheessa yksilöluetteloa laajennettiin lisäämällä siihen lampaiden ja vuohien maailmanlaajuisiin rotupaneeleihin (wpOA ja wpCH) kuuluvat eläimet sekä tässä vaiheessa saataville tulleet luonnonvaraiset lisänäytteet (4 O. orientalis ja 4 C. aegagrus). Genotyyppejä parannettiin ja vaiheistettiin ryhmien sisällä Beagle 453:lla, minkä jälkeen ne suodatettiin pois, jos genotyypin todennäköisyys oli alle 0,95. Lopuksi suodatimme pois kohteet, jotka olivat monomorfisia tässä tutkimuksessa käytettyjen yksilöiden eri osajoukoissa (ks. jäljempänä).

Vertaillaksemme Ovisin ja Capran välillä havaittuja valintasignaaleja suoritimme ristiintaulukoinnin kahden referenssigenomin välillä. Ensin käytimme Ensembl release 69 -koodipohjan54 pareittaista kohdistusputkea yhdenmukaistaaksemme lampaan (OARv3.1) ja vuohen (CHIR1.0) referenssigenomit. Tämä putki käyttää LastZ:tä55 DNA-tason kohdistamiseen, minkä jälkeen suoritetaan jälkikäsittely, jossa kohdistetut lohkot ketjutetaan toisiinsa niiden sijainnin mukaan molemmissa genomeissa. Ensembl ajaa LastZ-pareittaista kohdistusputkea rutiininomaisesti kaikille tuetuille lajeille, mutta vuohi ei ole vielä mukana Ensemblissä. Jotta vältettäisiin vääristymät kummankaan lajin suuntaan, tuotimme kaksi erilaista lajien välistä kohdistusta. Toisessa käytettiin lammasta referenssigenomina ja vuohta ei-referenssinä, kun taas toisessa käytettiin vuohta referenssigenomina ja lammasta ei-referenssinä. Ero on siinä, että referenssilajin genomialueiden on pakko karttua yksiselitteisesti ei-referenssilajin yksittäisiin lokuksiin, kun taas ei-referenssilajin genomialueiden sallitaan karttua useisiin referenssilajin paikkoihin. Saimme yhden referenssigenomin kromosomien segmenteille koordinaatit ei-referenssigenomissa. Lopuksi käytimme yhdessä suvussa havaittujen SNP:iden osalta koko genomin kohdistusta toisen suvun referenssigenomin kanssa vastaavien paikkojen tunnistamiseksi (lisätaulukko 6).

Genetic structure

Kuvaamaan ryhmien sisäistä geneettistä monimuotoisuutta käytimme VCFtools56 -ohjelmaa laskeaksemme geneettisen variaation yhteenvetotilastot 73 yksilölle Ovisin osalta (ts, 13 IROO, 20 IROA, 20 MOOA ja 20 wpOA) ja 72 yksilöä Capran osalta (eli 18 IRCA, 20 IRCH, 20 MOCH ja 14 wpCH). Mitatut tilastot olivat polymorfisten varianttien kokonaismäärä (S) koko yksilöjoukon osalta kussakin suvussa ja kussakin ryhmässä, keskimääräinen nukleotidididiversiteetti (π) kussakin ryhmässä ja sisäsiitoskerroin (F) kunkin yksilön osalta. Kunkin suvun sisällä luonnonvaraisen ryhmän ja kunkin kotieläinryhmän väliset erot testattiin käyttämällä yksipuolista t-testiä yksilöiden sisäsiitos- ja geneettisen kuorman arvojen osalta ja kaksipuolista Mann-Whitneyn testiä nukleotidididiversiteetin osalta paikkaa kohti.

Kunkin suvun sisällä neljän ryhmän välinen kokonaisdiversiteetti (ts, luonnonvaraiset, iranilaiset ja marokkolaiset kotieläimet sekä maailman paneeli) arvioitiin käyttämällä kaikkia biallelisia SNP:itä ja keskimääräistä painotettua pareittaista Fst:ää Weirin ja Cockerhamin57 mukaisesti, kuten VCFtoolsissa56 on toteutettu. Ryhmien välistä geneettistä rakennetta arvioitiin klusterointimenetelmällä sNMF26 sen jälkeen, kun aineistosta oli poistettu SNP:t, joiden yhteysepätasapaino (r²) oli yli 0,2 VCFtools-ohjelmalla. Linkage disequilibrium (r²) laskettiin SNP-parien välillä 50 SNP:n liukuvissa ikkunoissa, joista yksi SNP per pari poistettiin satunnaisesti, kun r² oli yli 0,2. Kutakin sNMF-analyysiä varten tehtiin viisi ajoa, joissa klusterien määrä (K) oli sama ja K:n arvot olivat 1-10. Käytimme risti-entropia-kriteeriä tunnistamaan todennäköisimmän klusterointiratkaisun, mutta tutkittiin myös vaihtoehtoisia jakoja eri K-luvuilla, jotta voitiin arvioida, miten yksilöt jakautuivat klustereiden välillä.

Jaetun syntyperän ja sekoittumisen erottamiseksi toisistaan ajettiin TreeMix-analyysi27 arvioimaan yhdessä populaatioiden jakautumisia ja myöhempiä sekoittumistapahtumia sNMF:ssä käytetyn karsitun aineiston avulla. Ajoimme TreeMix-ohjelman -global-optiolla tarkentaaksemme maksimaalisen todennäköisyyden johtopäätöksiä. Juurrutimme TreeMix-puun luonnonvaraisten ja kotieläinyksilöiden väliseen jakoon. Jackknifingin lohkokoko oli -k 500 SNP:tä, mikä vastaa suunnilleen 150 kb:tä, mikä ylittää sekä lampailla että vuohilla havaitut keskimääräiset LD-lohkot. Loimme Maximum Likelihood -puun, jossa ei ollut migraatiota, ja lisäsimme sitten migraatiotapahtumia ja tarkastelimme mallin selittämän varianssin ja yksilöiden välisten jäännösarvojen lisämuutosta. Tavoitteena oli havaita mahdolliset korkeat jäännösarvot tai vaellusreunat luonnonvaraisten ja kotieläinyksilöiden välillä. Tutkiaksemme tarkemmin TreeMix-ohjelmalla tunnistettujen mahdollisten sekoittumisvektoreiden tilastollista merkitystä (lisätaulukko 3) laskimme kolmen populaation testin f328 geneettisen introgression muodollisena testinä käyttäen ADMIXTOOLS-paketin qp3Pop-ohjelmaa58 kullekin ryhmäyhdistelmälle. Capran osalta wpCH-ryhmä jaettiin australialaisiin rotuihin, ranskalaisiin rotuihin ja italialaisiin rotuihin. Tulokset on raportoitu Supplementary Data 2:ssa.

Demografinen päättely

Kunkin suvun osalta suoritimme esivanhempien demografisen päättelyn analyysit käyttäen MSMC2-ohjelmistossa25 toteutettua MSMC-mallia. MSMC perustuu pareittain sekventiaaliseen markovilaiseen yhteensulautumiseen59; se käyttää kuitenkin syötteenä vaiheittaisten genomisekvenssitietojen haplotyyppejä. Kussakin analyysissä käytettiin kahta yksilöä yhdestä ryhmästä, joten haplotyyppejä oli neljä. Jokainen analyysi toistettiin toiselle satunnaiselle kahden yksilön joukolle, eli analyysi toistettiin ryhmäkohtaisesti. Syöttö- ja tulostiedostot luotiin ja analysoitiin MSMC-ohjelmiston mukana toimitetuilla python-skripteillä, jotka löytyvät osoitteesta https://github.com/stschiff/msmc-tools. Analyysiparametrit pidettiin oletusarvoisina lukuun ottamatta mutaationopeutta, joka asetettiin 2,5 × 10-8:aan, ja sukupolven pituutta, joka asetettiin 2 vuoteen. Arvioidaksemme aika-arvioiden epävarmuutta varioimme näitä parametreja (mutaationopeus 2,5×10-8 ja 1,0×10-8 yhdessä sukupolven pituuden ollessa 2 ja 4 vuotta) ja annoimme karkean arvion kesyyntymisajasta (ks. täydentävä kuva 2).

Geneettinen kuormitus

Geneettistä kuormitusta arvioitiin kahdella tavalla. Ensinnäkin laskemalla geneettinen kuormitus kullekin yksilölle kaikkien proteiineja koodaavien genomipositioiden haitallisten kuntovaikutusten summana Libradon ym.60 menetelmän mukaisesti. Lyhyesti sanottuna käytimme evolutiivisen rajoituksen korvikkeena nisäkkäiden 46-suuntaisesta kohdistuksesta saatuja PhyloP-pisteitä (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). Tästä linjauksesta tunnistimme proteiineja koodaavat paikat, jotka kehittyvät toiminnallisten rajoitusten alaisina (phyloP-pisteet ≥1,5). Kunkin Ovis- tai Capra-genomin osalta tutkittiin sitten, olivatko nämä kohdat mutatoituneet. Jos näin oli, laskimme phyloP-pisteet yhteen kaikkien mutaatiokohtien osalta niin, että mutaatiot erittäin rajoitetuissa kohteissa vaikuttavat suhteellisesti enemmän kokonaiskuormitusarvioon. Näin saatiin kuormitusarvio kullekin lampaan/vuohen genomille. Lopuksi keskimääräisen kuormituksen saamiseksi sivustokohtaisesti jaoimme sen analysoitujen paikkojen kokonaismäärällä. On syytä huomata, että ehdollistimme homotsygoottiset paikat välttyäksemme heterotsygoottisten paikkojen mutaatioiden dominanssikertoimen mallintamiselta (esim. resessiivinen, intermediäärinen, dominoiva). Toiseksi vertasimme geenikohtaisesti geneettistä haitallista kuormitusta luonnonvaraisissa ja kotieläiminä pidetyissä Ovis-ryhmissä suorittamalla Wilcoxonin testin, jossa vaihtoehtoinen hypoteesi oli, että kotieläimillä on enemmän kuormitusta kuin luonnonvaraisilla sukulaisilla. p-arvot korjattiin moninkertaisen testauksen vuoksi61 , ja käytimme kynnysarvoa, jonka mukaan oikaistujen p-arvojen on oltava < 0,05. Suoritimme geeniontologian rikastusanalyysin niille geeneille, jotka osoittivat merkittävää geneettisen kuormituksen lisääntymistä WebGestalt-ohjelmalla62,63. Koska referenssigenomit ovat huonosti annotoituja geenien osalta, tukeuduimme yhden kopion ortologeihin lajimme sekä ihmisen ja hiiren välillä. X-kromosomin geenit jätettiin taustajoukon ulkopuolelle. Emme tehneet tätä analyysia Capralle, koska luonnonvaraisissa näytteissä havaittiin suurempaa sisäsiittoisuutta.

Valinnan tunnusmerkkien havaitseminen

Kotieläinjalostukseen liittyvän valinnan tunnusmerkkien havaitsemiseksi käytimme kaikkia biallelisia SNP:itä, joiden pieni alleelifrekvenssi oli suurempi kuin 0,10 vähintään yhdessä kolmesta testatusta ryhmästä (eli iranilaisista ja marokkolaisista kotieläinjalostajaryhmistä sekä kunkin suvun luonnonvaraisesta ryhmästä). Koska odotimme, että domestikaatioprosessiin liittyvän valinnan merkkejä esiintyisi kaikissa kotieläimissä, noudatimme seuraavaa yleistä strategiaa: testasimme hapFLK29:llä (ks. lisähuomautus 5 ja lisäkuviot 9, 10 ja 11) kunkin suvun osalta luonnonvaraista ryhmää kutakin perinteisesti hoidettua kotieläinryhmää vastaan (eli iranilaista ja marokkolaista ryhmää) ja keskityimme niihin yhteisiin alueisiin, jotka olivat oletettavasti valinnan kohteena ja jotka havaittiin molemmissa tapauksissa. Ryhmien otoskoko (n = 13-20) oli menetelmän vaatimusten mukainen29. Tarkistimme silmämääräisesti, oliko hapFLK:n avulla havaittuja johdonmukaisia valinnan merkkejä myös kunkin suvun vastaavassa maailmanlaajuisessa paneelisarjassa, mutta emme ottaneet näitä ryhmiä mukaan tilastolliseen testiin niiden monirotuisen koostumuksen vuoksi. Lopuksi etsimme yhteisiä valinnan merkkejä Ovisin ja Capran välillä käyttämällä ositettua FDR-menetelmää. Strategia on esitetty täydentävässä kuvassa 4.

Toteutimme hapFLK-testejä luonnonvaraisen ryhmän ja kunkin suvun kunkin iranilaisen ja marokkolaisen ryhmän erottamiseksi toisistaan. Sukulaisuusmatriisi laskettiin ryhmäparien välisistä Reynoldin geneettisistä etäisyyksistä64 käyttäen satunnaista yhden prosentin osajoukkoa varianteista. Johdettu populaatiopuu rakennettiin käyttämällä naapuriliitosalgoritmia. Jokaiselle SNP:lle suoritettiin hapFLK-testi, joka sisältää haplotyyppistä tietoa valikoivien pyyhkäisyjen havaitsemisen tehon lisäämiseksi. Kunkin testatun SNP:n osalta hapFLK-statistiikka laski haplotyyppisten frekvenssien poikkeaman sukulaisuusmatriisin avulla estimoituun neutraaliin malliin nähden65. Linkage disequilibrium -tiedon hyödyntämiseksi hapFLK käyttää Scheetin ja Stephensin66 monipistemallia monilokusgenotyypeille, joka voidaan sovittaa vaiheistamattomaan dataan. Yksi tämän mallin tärkeimmistä sovelluksista on vaiheen estimointi (fastPHASE-ohjelmisto66). Analyysissämme malli koulutettiin vaiheistamattomilla tiedoilla, joten analyysissämme otetaan huomioon vaiheeseen liittyvä epävarmuus. Menetelmää käytettiin paikallisten haplotyyppien ryhmittelyyn kromosomeja pitkin määritettyyn klusterien määrään K, joka asetettiin 25:ksi, käyttäen Piilotetun Markovin mallia.

Tunnistaaksemme yhteiset alueet, jotka oletettavasti ovat valinnan kohteena kahdessa perinteisesti hoidetussa kotieläinryhmässä kunkin suvun osalta, yhdistimme kaksi aiempaa hapFLK-analyysiä. Kunkin analyysin osalta hapFLK-pisteet sovitettiin χ2-jakaumaan p-arvojen saamiseksi (käsikirjoitus saatavilla osoitteessa https://forge-dga.jouy.inra.fr/projects/hapflk/documents). Kahden luonnonvaraisen ryhmän ja kunkin kotieläinryhmän välisen kontrastin tulokset yhdistettiin Stoufferin menetelmällä67 , jotta saatiin yksittäiset p-arvot luonnonvaraisten ja kotieläinten väliselle vertailulle. Lopuksi koko SNP-joukkoon sovellettiin FDR-menetelmää68 yhdistettyjen p-arvojen muuntamiseksi q-arvoiksi. SNP:t, joiden q-arvot olivat < 10-2, säilytettiin ja ryhmiteltiin genomialueiksi, kun ne olivat alle 50 kb:n etäisyydellä toisistaan.

Tullaksemme selville, oliko valinnan signaali yhteinen Ovisin ja Capran välillä, käytimme aluksi kahden referenssigenomin ristiinkohdistusta homologisten segmenttien tunnistamiseksi. Sen jälkeen sovelsimme ositettua FDR-kehystä69. Tämä lähestymistapa perustuu siihen, että testeissä on luonnostaan kerrostuneisuutta, kun otetaan huomioon geneettisen datan ennakkotiedot69 , koska todellisten vaihtoehtoisten hypoteesien taustalla oleva jakauma voi olla erilainen eri genomialueiden erilaisen dynamiikan mukaan, mikä johtaa erilaisiin p-arvojen jakaumiin. Tämä edellyttää FDR-korjattujen p-arvojen (eli q-arvojen) saamista erikseen eri ositteille. Etsimme yhteneväisyyksiä kussakin suvussa erottelemalla toisesta suvusta havaittujen alueiden kanssa homologiset alueet (jota kutsutaan yhteiseksi kerrostumaksi) ja loput genomista (jota kutsutaan yleiseksi kerrostumaksi). Otimme p-arvot erikseen kummallekin määritellylle kerrokselle ja laskimme sitten q-arvot FDR-menetelmän avulla. Nämä stratifioidut q-arvot olivat lopullisia suureita, joita tarkasteltiin tilastollisen merkitsevyyden (<10-2) kannalta, jotta valinnan kohteena olevat SNP:t voitiin havaita ja yhdistää vastaaviin genomialueisiin.

Testataksemme luonnonvaraiset eläimet ja kotieläimet erottavan valinnan yhteneviä merkkejä molemmissa suvuissa, tarkastelimme q-arvoihin sovelletun merkitsevyyskynnyksen (jonka teimme vaihtelevaksi 0,2:sta 0,002:een) välistä suhdetta jommassakummassa suvussa ja arvioitua todennäköisyyttä sille, että SNP on valikoitunut toisen suvun yhteisessä stratumissa, käyttäen Storeyn ym. menetelmää.70 lähestymistapa. Laskennallisen todennäköisyyden kasvu q-arvoon sovelletun kynnysarvon pienentyessä (tiukuuden lisääntyessä) osoittaa, että mitä merkittävämpi alue on yhdessä suvussa, sitä todennäköisemmin löydämme merkitseviä SNP:itä toisesta suvusta.

Suodatimme pois valikoitumissignaalit, jotka eivät olleet johdonmukaisia kolmen kotimaisten ryhmien välillä. Käytimme kunkin havaitun alueen osalta yksilöiden vaiheistettuja haplotyyppejä, jotka klusteroitiin käyttämällä Neighbor-Joining-puita sekvenssien välisen identiteettiprosentin perusteella. Ainoastaan alueet, jotka osoittivat johdonmukaisia signaaleja, säilytettiin (Täydentävä kuva 5).

Johtopäätökseksi siitä, viittasivatko hapFLK:lla havaitut valintasignaalit valinnan relaksaatioon vai positiiviseen valintaan kotieläimissä, arvioimme eron nukleotididiversiteetissä (π) kullakin oletetulla valinnan kohteena olevalla alueella luonnonvaraisten ja kotieläinryhmien välillä. Ilmaisimme tämän eron Δπ-indeksinä, joka laskettiin kullekin genomialueelle luonnonvaraiselle ryhmälle lasketun π:n ja iranilaisille ja marokkolaisille kotieläinryhmille lasketun π:n keskiarvon erotuksena vähennettynä näiden kahden ryhmän välisellä π:n erotuksella, joka on laskettu koko genomille:

Negatiivinen arvo osoittaisi, että nukleotidien monimuotoisuus on alhaisempi villiryhmässä verrattuna kahden kotieläinryhmän keskiarvoon, ja sen katsottaisiin osoittavan valinnan relaksoitumista näissä viimeksi mainituissa ryhmissä, monipuolistavaa valintaa kotieläimissä tai positiivista valintaa villeissä. Sitä vastoin positiivinen arvo osoittaisi kotieläinryhmissä tapahtunutta suuntautunutta positiivista tai vakauttavaa valintaa. Haplotyyppiklusteroinnin avulla tarkistimme myös manuaalisesti kullakin alueella, oliko havaittu valikoiva pyyhkäisy vahvistanut Δπ-indeksin antamat viitteet.

Suoritimme toiminnallisia tulkintoja seuraavasti. Kunkin valinnan kohteena olevan alueen osalta tarkastelimme aluetta plus 50 kb kummallakin puolella toiminnallisten roolien tunnistamiseksi ja 5 kb ylävirtaan ja alavirtaan geeneistä ja arvioimme näiden koordinaattien päällekkäisyyttä kiinnostavien geenien säilyttämiseksi. Lopuksi katsoimme, että geeni liittyi tiettyyn havaittuun alueeseen, kun alueen ja geenin sijainnit olivat päällekkäisiä. Tämän jälkeen arvioimme, mikä geeni oli todennäköisimmin valinnan kohteena, tarkastelemalla huippusignaalia lähimpänä olevaa geeniä eli pienimmän q-arvon sijaintia alueen sisällä. Geenit annotoitiin toiminnallisesti Uniprotin (http://www.uniprot.org/) avulla ottamalla huomioon niiden osallistuminen ”Biologinen prosessi” -luokan (eli GO:0008150) 30 lapsitermiin (eli termien suoriin aleneviin). Haimme kaikki kutakin geeniä vastaavat GO-termit (Supplementary Data 4) 30:lle 33 kategoriasta, koska emme ottaneet huomioon kolmea termiä, jotka eivät osallistuneet nisäkkäiden toimintoihin (ts. GO:0006791 rikin hyödyntäminen, GO:0006794 fosforin hyödyntäminen, GO:0015976 hiilen hyödyntäminen). Suoritimme kaksi χ2-testiä vertaillaksemme geenien jakaumia GO-luokissa, eli (i) lajikohtaisilta alueilta valikoituneita geenejä verrattuna homologisilta alueilta valikoituneisiin geeneihin ja (ii) kaikkia valikoituneita geenejä verrattuna 18 689 ihmisen geeneihin, jotka oli liitetty GO-termeihin Swiss-Protissa. Tulkitaksemme geenien toimintoja kotieläinkontekstissa haimme myös kirjallisuudesta saatavilla olevaa tietoa niiden fenotyyppisistä vaikutuksista.

Viimeiseksi löysimme aiemmin havaituilta alueilta ne SNP:t, jotka erosivat eniten luonnonvaraisten ja kotieläinryhmien välillä, käyttämällä FLK-statistiikkaa. Kuten hapFLK, se edustaa yksittäisten merkkien alleelifrekvenssien poikkeamaa suhteessa sukulaisuusmatriisin avulla estimoituun neutraaliin malliin65. Kahden analyysin tulokset sovitettiin χ2-jakaumaan ja saadut p-arvot yhdistettiin samalla menettelyllä kuin hapFLK-testissä. P-arvojen epätasainen jakauma esti kuitenkin FDR-kehyksen soveltamisen, ja valitsimme hapFLK:lla havaituilla alueilla olevat SNP:t, joiden p-arvot olivat <10-4. Näiden SNP:iden osalta käytimme Variant Effect Predictor (VEP) -annotaatioita71 , jotka luotiin Ensembl v74 -lampaan OARv3.1 -genomin annotaatiosta Ovisille (http://www.ensembl.org/Ovis_aries/Tools/VEP) ja NCBI:n eukaryoottisen genomin annotaatioputken tuottamasta vuohen CHIR1.0 -genomin annotaatiosta Capralle (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/). SNP:t luokiteltiin intergeenisiin, ylävirtaan ja alavirtaan (mukaan lukien UTR:t) sekä intronisiin ja eksonisiin asemiin. Erot niiden SNP:iden jakaumien välillä, joiden FLK:n p-arvo oli <10-4, ja kaikkien valikoitumissignatuurien havaitsemiseen käytettyjen SNP:iden välillä tutkittiin χ2-testillä.

Tietojen saatavuus

Sekvenssi- ja metatietoaineistot, jotka on tuotettu 73:sta näissä analyyseissä käytetystä Ovis-näytteestä ja 72:sta Capra-näytteestä, ovat julkisesti saatavilla. Yleiset tiedot ja kaikki vcf-tiedostot löytyvät Ensemblin verkkosivuilta (http://projects.ensembl.org/nextgen/). Kaikki Fastq-tiedostot, Bam-tiedostot ja O. orientaliksen ja C. aegagrusin de novo -assemblaatiot löytyvät European Nucleotide Archive -sivustolta (https://www.ebi.ac.uk/ena) Nextgen-projektin liittymiskoodilla (PRJEB7436).