Campionamento

Pecore domestiche (O. aries) e capre (C. hircus) sono state campionate in Iran (gruppi IROA e IRCH, rispettivamente) e Marocco (gruppi MOOA e MOCH, rispettivamente) per un totale di 20 animali per gruppo (Fig. 6 supplementare). Questi campioni sono stati raccolti tra gennaio 2008 e marzo 2012 nella parte settentrionale del Marocco e tra agosto 2011 e luglio 2012 nell’Iran nord-occidentale, nell’ambito del progetto europeo Nextgen (Grant Agreement no. 244356) in conformità alle norme etiche della direttiva 86/609/CEE dell’Unione europea. I fermagli auricolari sono stati raccolti dalla parte distale dell’orecchio di animali scelti a caso, e immediatamente conservati in etanolo al 96% per un giorno prima di essere trasferiti in perle di gel di silice fino all’estrazione del DNA.

Le specie selvatiche muflone asiatico (O. orientalis) e stambecco Bezoar (C. aegagrus) sono state campionate nell’Iran nord-occidentale all’interno della culla dell’addomesticamento21,22. Tredici mufloni asiatici e 18 tessuti di stambecco Bezoar (rispettivamente, gruppi IROO e IRCA, Fig. 6 supplementare) sono stati raccolti da animali in cattività o cacciati di recente, e da campioni congelati disponibili presso il Dipartimento dell’ambiente iraniano. Questo approccio di campionamento basato sull’individuo è stato progettato per minimizzare i potenziali errori evitando la sovrarappresentazione degli effetti locali (ad es,

Dati aggiuntivi

Inoltre, è stato assemblato un pannello di razza mondiale per pecore e capre (wpOA e wpCH, rispettivamente). wpOA comprendeva 20 campioni di ri-sequenziamento del genoma intero (WGS) a copertura 12x che rappresentano 20 diverse razze mondiali fornite dall’International Sheep Genome Consortium. wpCH consisteva di 14 campioni WGS sequenziati a copertura 12x che rappresentano 9 individui europei, cioè, 2 campioni francesi alpini e 2 campioni francesi Saanen sequenziati dall’INRA, 5 campioni italiani Saanen forniti dal Parco Tecnologico Padano, e 5 individui australiani, cioè, 2 Boer, 2 Rangeland e 1 Cashmere forniti dal CSIRO (Dati supplementari 5).

Produzione di dati WGS

Il DNA genomico è stato estratto con successo da tutti i campioni di tessuto utilizzando il kit Macherey Nagel NucleoSpin 96 Tissue, adattando il protocollo del produttore. Il campionamento del tessuto è stato eseguito in blocchi di pozzetti quadrati MN per ottenere 25 mg di frammenti per campione. Tre e mezzo MN piazza-96 blocchi sono stati preparati, e l’estrazione è stata eseguita utilizzando un Tecan Freedom Evo Liquid handler seguendo il protocollo del produttore. Una fase di pre-lisi è stata effettuata per omogeneizzare i campioni con 180 µl di T1 Buffer e 25 µl di proteinasi K per una notte a 56 °C. Per regolare le condizioni di legame, sono stati aggiunti 200 µl di buffer BQ1 e la piastra dei campioni è stata incubata 1 ora a 70 °C; successivamente sono stati aggiunti 200 µl di etanolo al 100%. I lisati sono stati trasferiti alla piastra di legame Nucleospin Tissue e un vuoto (-0,2 bar, 5 min) è stato applicato per rimuovere il flusso passante. Tre passaggi di lavaggio sono stati fatti con BW e B5 buffer, rispettivamente, e un vuoto è stato applicato di nuovo per scartare il flusso attraverso. Prima dell’eluizione del DNA genomico, una membrana di silice Nucleospin Tissue binding plate è stata asciugata sotto vuoto con almeno -0,6 bar per 10 minuti. La fase di eluizione è stata eseguita con 100 µl di tampone BE preriscaldato (70 °C) e una centrifugazione a 3700 rcf per 5 min nei pozzetti 96-PCR. Il DNA genomico è stato conservato a 4 °C per evitare il congelamento e testato per la concentrazione (come ng/μl) utilizzando il metodo Picogreen e utilizzando un Nanodrop.

I genomi interi sono stati risequenziati da 500 ng di DNA genomico che sono stati tosati a un intervallo di 150-700 bp utilizzando lo strumento Covaris® E210 per ogni campione e utilizzati per la preparazione della libreria Illumina® con un protocollo semi-automatizzato. Riparazione finale, A-tailing, e Illumina ® adattatori compatibili (BioScientific) legatura sono stati eseguiti utilizzando il SPRIWorks Library Preparation System e SPRI TE strumento (Beckmann Coulter) seguendo il protocollo del produttore. Una selezione delle dimensioni di 300-600 bp è stata applicata per recuperare la maggior parte dei frammenti. I frammenti di DNA sono stati amplificati da 12 cicli di PCR utilizzando Platinum Pfx Taq Polymerase Kit (Life® Technologies) e Illumina® adapter-specific primers. Le librerie sono state purificate con perline AMPure XP 0.8x (Beckmann Coulter), e analizzate con il Bioanalyzer Agilent 2100 (Agilent® Technologies) e la quantificazione qPCR. Le librerie sono state sequenziate utilizzando la chimica di lettura a 100 basi in cella di flusso paired-end su Illumina® HiSeq2000.

Le letture paired-end di Illumina per Ovis sono state mappate sul genoma di riferimento della pecora (OAR v3.1, assemblaggio GenBank GCA_000298735.146), e per Capra sul genoma di riferimento della capra (CHIR v1.0, assemblaggio GenBank GCA_000317765.147), utilizzando BWA-MEM48. Il file BAM prodotto per ogni individuo è stato ordinato utilizzando Picard SortSam e migliorato utilizzando sequenzialmente Picard MarkDuplicates (http://picard.sourceforge.net), GATK RealignerTargetCreator e GATK IndelRealigner49, e SAMtools calmd50.

La scoperta delle varianti è stata effettuata utilizzando tre diversi algoritmi: Samtools mpileup50, GATK UnifiedGenotyper51, e Freebayes52. I siti delle varianti sono stati identificati indipendentemente per ciascuno dei sei gruppi, utilizzando le modalità multi-campione degli algoritmi di chiamata: (i) 162 campioni da MOOA; (ii) 20 campioni da IROA; (iii) 14 campioni da IROO; (iv) 162 campioni da MOCH; (v) 20 campioni da IRCH; (vi) 19 campioni da IRCA. Per alcuni gruppi, le WGS di più individui erano disponibili come parte del progetto NextGen (vedi sopra). I campioni utilizzati nel presente studio sono stati selezionati per ottenere gruppi bilanciati di 20 individui, quando possibile. Per i gruppi IRCA e IROO, ulteriori campioni sono diventati disponibili in una fase successiva e sono stati aggiunti per le analisi a valle. Gli animali con bassa qualità di allineamento e di chiamata sono stati rimossi per ottenere il set di dati finale (dati supplementari 5).

In ogni gruppo, ci sono stati due giri successivi di filtraggio della qualità del sito delle varianti. La fase di filtraggio 1 ha unito le chiamate dei tre algoritmi, filtrando le chiamate a più bassa confidenza. Un sito variante è passato se è stato chiamato da almeno due diversi algoritmi di chiamata con qualità della variante phred >30. Un allele alternativo in un sito è passato se è stato chiamato da uno qualsiasi degli algoritmi di chiamata, e il conteggio del genotipo era >0. La fase di filtraggio 2 ha utilizzato la ricalibrazione del punteggio di qualità delle varianti di GATK. In primo luogo, abbiamo generato un set di allenamento dei siti di varianti a più alta confidenza all’interno del gruppo in cui (i) il sito è chiamato da tutti e tre i chiamanti di varianti con phred-scaled variant quality >100, (ii) il sito è biallelico, (iii) il conteggio degli alleli minori è almeno 3, contando solo i campioni con phred-scaled genotype quality >30. Il set di allenamento è stato utilizzato per costruire un modello gaussiano utilizzando lo strumento GATK VariantRecalibrator utilizzando le seguenti annotazioni di varianti da UnifiedGenotyper: QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. Un modello gaussiano è stato applicato all’intero set di dati, generando un VQSLOD (log odds ratio di essere una vera variante). I siti sono stati filtrati se VQSLOD <valore di cutoff. Il valore di cutoff è stato impostato per ogni gruppo come segue: VQSLOD minimo = {il valore mediano di VQSLOD per le varianti del training set}-3 × {la deviazione assoluta mediana VQSLOD delle varianti del training set}. Il rapporto SNP di transizione/trasversione ha suggerito che il criterio di cutoff scelto ha dato il miglior equilibrio tra selettività e sensibilità.

Sono stati generati set di chiamate SNPs per sei gruppi di animali Ovis e Capra (cioè, domestici iraniani e marocchini, e selvatici per ogni genere). Poiché le analisi eseguite in questo studio richiedevano confronti intergruppo, abbiamo creato set di chiamate di genotipi in un insieme coerente di siti SNP per tutti gli animali di ogni gruppo. Per ogni genere, abbiamo unito i siti di chiamata delle varianti dei suoi tre gruppi e abbiamo mantenuto solo le posizioni bialleliche senza dati mancanti. I genotipi sono stati ri-chiamati in ogni sito SNP biallelico per tutti gli individui di interesse da GATK UnifiedGenotyper utilizzando l’opzione GENOTYPE_GIVEN_ALLELES. In questa fase, l’elenco degli individui è stato ampliato per includere gli animali appartenenti ai panel mondiali delle razze ovine e caprine (wpOA e wpCH) e ulteriori campioni selvatici che si sono resi disponibili in questa fase (4 O. orientalis e 4 C. aegagrus). I genotipi sono stati migliorati e sfasati all’interno dei gruppi da Beagle 453, e poi filtrati dove la probabilità del genotipo era inferiore a 0,95. Infine, abbiamo filtrato i siti che erano monomorfi attraverso i diversi sottoinsiemi di individui utilizzati in questo studio (vedi sotto).

Al fine di confrontare i segnali di selezione rilevati tra Ovis e Capra, abbiamo eseguito un allineamento incrociato tra i due genomi di riferimento. In primo luogo, abbiamo usato la pipeline di allineamento a coppie dal codice base Ensembl release 6954 per allineare i genomi di riferimento di pecora (OARv3.1) e capra (CHIR1.0). Questa pipeline utilizza LastZ55 per allineare a livello di DNA, seguita da una post-elaborazione in cui i blocchi allineati sono concatenati in base alla loro posizione in entrambi i genomi. La pipeline di allineamento a coppie LastZ viene eseguita di routine da Ensembl per tutte le specie supportate, ma la capra non è ancora inclusa in Ensembl. Per evitare distorsioni verso una delle due specie, abbiamo prodotto due diversi allineamenti interspecifici. Uno ha usato la pecora come genoma di riferimento e la capra come non riferimento, mentre l’altro ha usato la capra come genoma di riferimento e la pecora come non riferimento. La differenza è che le regioni genomiche della specie di riferimento sono costrette a mappare in modo univoco su singoli loci della specie non di riferimento, mentre le regioni genomiche non di riferimento possono mappare su più posizioni della specie di riferimento. Abbiamo ottenuto per i segmenti di cromosomi di un genoma di riferimento le coordinate sul genoma non di riferimento. Infine, per gli SNPs scoperti in un genere, abbiamo utilizzato l’allineamento dell’intero genoma con il genoma di riferimento dell’altro genere per identificare le posizioni corrispondenti (Tabella 6 supplementare).

Struttura genetica

Per descrivere la diversità genetica all’interno dei gruppi, abbiamo utilizzato VCFtools56 per calcolare le statistiche riassuntive della variazione genetica sui 73 individui per Ovis (cioè, 13 IROO, 20 IROA, 20 MOOA, e 20 wpOA) e 72 individui per Capra (cioè 18 IRCA, 20 IRCH, 20 MOCH, e 14 wpCH). Le statistiche misurate erano il numero totale di varianti polimorfiche (S) per l’intero set di individui in ogni genere e all’interno di ogni gruppo, la diversità nucleotidica media (π) all’interno di ogni gruppo e il coefficiente di inbreeding (F) per ogni individuo. All’interno di ogni genere, le differenze tra il gruppo selvatico e ogni gruppo domestico sono state testate utilizzando un test t unilaterale per i valori individuali di inbreeding e carico genetico, e un test Mann-Whitney bilaterale per la diversità nucleotidica per sito.

La divergenza complessiva tra i quattro gruppi all’interno di ogni genere (cioè, selvatici, domestici iraniani e marocchini e panel mondiale) è stata stimata utilizzando tutti gli SNP biallelici e la media ponderata di Fst a coppie seguendo Weir e Cockerham57 come implementato in VCFtools56. La struttura genetica tra i gruppi è stata valutata con il metodo di clustering sNMF26, dopo la potatura del set di dati per rimuovere SNPs con linkage disequilibrium (r²) maggiore di 0,2 utilizzando VCFtools. Linkage disequilibrium (r²) è stato calcolato tra coppie di SNPs all’interno di finestre scorrevoli di 50 SNPs, con uno SNP per coppia rimosso a caso quando r² era maggiore di 0.2. Per ogni analisi sNMF, sono state eseguite cinque corse dello stesso numero di cluster (K) con valori di K da 1 a 10. Abbiamo usato il criterio dell’entropia incrociata per identificare la soluzione di clustering più probabile, tuttavia, sono state esplorate anche partizioni alternative per diversi numeri di K per valutare come gli individui sono stati divisi tra i cluster.

Per distinguere tra ascendenza condivisa e commistione, abbiamo eseguito TreeMix27 per stimare congiuntamente le suddivisioni della popolazione e i successivi eventi di commistione utilizzando il set di dati sfrondato usato per sNMF. Abbiamo eseguito TreeMix con l’opzione -global per raffinare le nostre inferenze di massima verosimiglianza. Abbiamo radicato l’albero TreeMix con la divisione tra individui selvatici e domestici. La dimensione del blocco per il jackknifing era -k 500 SNPs, che corrisponde approssimativamente a 150 kb, superando i blocchi medi di LD trovati sia nelle pecore che nelle capre. Abbiamo generato un albero di Maximum Likelihood senza migrazione e poi abbiamo aggiunto eventi di migrazione ed esaminato il cambiamento incrementale nella varianza spiegata dal modello e i valori residui tra gli individui. L’obiettivo era quello di rilevare qualsiasi potenziale alto valore residuo o bordo di migrazione tra gli individui selvatici e domestici. Per esplorare ulteriormente la rilevanza statistica dei possibili vettori di commistione identificati da TreeMix (Tabella supplementare 3), abbiamo calcolato il test a tre popolazioni f328 come test formale di introgressione genetica, utilizzando il programma qp3Pop della suite ADMIXTOOLS58 per ogni combinazione di gruppi. Per Capra, il gruppo wpCH è stato diviso tra razze australiane, razze francesi e razze italiane. I risultati sono riportati nei dati supplementari 2.

Inferenza demografica

Per ogni genere, abbiamo effettuato analisi di inferenza demografica ancestrale utilizzando il modello MSMC implementato nel software MSMC225. MSMC si basa sulla coalescenza sequenziale markoviana a coppie59; tuttavia, utilizza come input gli aplotipi dei dati di sequenza del genoma a fasi. Per ogni analisi abbiamo usato due individui da un gruppo, quindi 4 aplotipi. Ogni analisi è stata ripetuta per un altro set casuale di due individui, cioè una replica dell’analisi per gruppo. I file di input e di output sono stati generati e analizzati con gli script python forniti con il software MSMC e che si trovano a https://github.com/stschiff/msmc-tools. I parametri di analisi sono stati mantenuti come predefiniti, tranne il tasso di mutazione che è stato impostato a 2,5×10-8 e la lunghezza della generazione è stata impostata a 2 anni. Al fine di stimare l’incertezza sulle stime temporali, abbiamo variato questi parametri (tasso di mutazione di 2,5×10-8 e 1,0×10-8 in combinazione con la lunghezza della generazione di 2 e 4 anni) e abbiamo fornito una stima approssimativa del periodo di domesticazione (vedi Fig. 2 supplementare).

Carico genetico

Il carico genetico è stato stimato in due modi. In primo luogo, calcolando il carico genetico per ogni individuo come la somma degli effetti deleteri di fitness su tutte le posizioni genomiche codificanti le proteine seguendo il metodo di Librado et al.60. Brevemente, come proxy per il vincolo evolutivo, abbiamo usato i punteggi PhyloP dall’allineamento dei mammiferi a 46 vie (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/). Da questo allineamento, abbiamo identificato i siti che codificano le proteine che evolvono sotto vincoli funzionali (punteggio PhyloP ≥1.5). Per ogni genoma di Ovis o Capra, abbiamo poi indagato se questi siti erano mutati. In tal caso, abbiamo sommato i punteggi phyloP su tutti i siti mutati, in modo che le mutazioni in siti altamente vincolati contribuiscano proporzionalmente di più alla stima del carico totale. Questo ha fornito una stima del carico per ogni genoma di pecora/capra. Infine, per ottenere un carico medio per sito, abbiamo diviso per il numero totale di posizioni analizzate. Vale la pena notare che abbiamo condizionato i siti omozigoti per evitare di modellare il coefficiente di dominanza delle mutazioni nei siti eterozigoti (ad esempio, recessivi, intermedi, dominanti). In secondo luogo, abbiamo confrontato gene per gene il carico genetico deleterio nei gruppi di Ovis selvatici e addomesticati eseguendo un test di Wilcoxon, con l’ipotesi alternativa che gli animali domestici hanno più carico dei parenti selvatici. I valori di p sono stati corretti per test multipli61 e abbiamo applicato una soglia di p-valori aggiustati < 0,05. Abbiamo eseguito un’analisi di arricchimento dell’ontologia genica sull’insieme dei geni che mostrano un aumento significativo del carico genetico usando WebGestalt62,63. Poiché i genomi di riferimento sono scarsamente annotati per i geni, ci siamo basati su ortologhi a copia singola tra la nostra specie e umano e mouse. I geni del cromosoma X sono stati esclusi dal set di sfondo. Non abbiamo effettuato questa analisi su Capra a causa della maggiore consanguineità osservata nei campioni selvatici.

Rilevamento delle firme di selezione

Per rilevare le firme di selezione relative all’addomesticamento, abbiamo utilizzato tutti gli SNP biallelici che mostravano una frequenza allelica minore maggiore di 0,10 in almeno uno dei tre gruppi testati (cioè, i gruppi domestici iraniani e marocchini, e il gruppo selvatico per ogni genere). Poiché ci aspettavamo che le firme della selezione legata al processo di addomesticamento fossero presenti in tutti gli animali domestici, abbiamo adottato la seguente strategia generale: abbiamo testato con hapFLK29 (vedi nota supplementare 5 e figure supplementari 9, 10 e 11) per ogni genere il gruppo selvatico contro ciascuno dei gruppi domestici gestiti tradizionalmente (cioè, l’iraniano e il Marocco) e ci siamo concentrati su quelle regioni comuni putativamente sotto selezione che sono state rilevate in entrambi i casi. Le dimensioni del campione del gruppo (n = 13-20) erano compatibili con i requisiti del metodo29. Abbiamo controllato visivamente se le firme coerenti di selezione trovate con hapFLK fossero presenti anche nel corrispondente set di pannelli mondiali di ogni genere, ma non abbiamo incluso questi gruppi nel test statistico a causa della loro composizione multirazza. Infine, abbiamo cercato segnali condivisi di selezione tra Ovis e Capra utilizzando un approccio FDR stratificato. La strategia è illustrata nella Fig. 4.

Abbiamo eseguito i test hapFLK per contrapporre il gruppo selvatico a ciascuno dei gruppi iraniani e marocchini in ogni genere. La matrice di parentela è stata calcolata dalle distanze genetiche di Reynold64 tra le coppie di gruppi, utilizzando un sottoinsieme casuale dell’uno per cento delle varianti. L’albero di popolazione dedotto è stato costruito usando l’algoritmo neighbor-joining. Per ogni SNP, abbiamo eseguito il test hapFLK che incorpora le informazioni aplotipiche per aumentare il potere di rilevare gli sweep selettivi. Per ogni SNP testato, la statistica hapFLK ha calcolato la deviazione delle frequenze aplotipiche rispetto al modello neutro stimato dalla matrice di parentela65. Per sfruttare le informazioni sul linkage disequilibrium, hapFLK utilizza il modello multipunto di Scheet e Stephens’66 per i genotipi multilocus che può essere adattato a dati non fasati. Una delle principali applicazioni di questo modello è quella di eseguire la stima della fase (software fastPHASE66). Nella nostra analisi, il modello è stato addestrato su dati non sfasati, e quindi la nostra analisi tiene conto dell’incertezza di fase. Il metodo è stato utilizzato per raggruppare gli aplotipi locali lungo i cromosomi in un numero specificato di cluster K impostato su 25, utilizzando un modello Hidden Markov Model.

Per identificare le regioni comuni putativamente sotto selezione nei due gruppi domestici tradizionalmente gestiti per ogni genere, abbiamo combinato le due analisi hapFLK precedenti. Per ogni analisi i punteggi hapFLK sono stati adattati a una distribuzione χ2 per ottenere i valori di p (script disponibile a https://forge-dga.jouy.inra.fr/projects/hapflk/documents). I risultati dei due contrasti tra il gruppo selvatico e ciascuno dei gruppi domestici sono stati combinati utilizzando il metodo di Stouffer67 per ottenere singoli p-valori per il confronto tra animali selvatici e domestici. Infine, la struttura FDR68 è stata applicata all’intero set di SNPs per convertire i valori p combinati in valori q. Gli SNP che mostravano valori q < 10-2 sono stati conservati e raggruppati in regioni genomiche quando erano distanti meno di 50 kb l’uno dall’altro.

Per indagare se il segnale di selezione era condiviso tra Ovis e Capra, abbiamo prima utilizzato l’allineamento incrociato dei due genomi di riferimento per identificare segmenti omologhi. Abbiamo poi applicato un quadro FDR stratificato69. Questo approccio si basa sul fatto che c’è una stratificazione intrinseca nei test date le informazioni precedenti nei dati genetici69, perché la distribuzione sottostante di vere ipotesi alternative potrebbe essere diversa a seconda delle diverse dinamiche di varie regioni genomiche, portando a diverse distribuzioni di p-valori. Questo richiede di ottenere p-valori aggiustati FDR (cioè, q-valori) separatamente per i diversi strati. Abbiamo cercato le convergenze in ogni genere separando le regioni omologhe a quelle rilevate nell’altro genere (indicato come strato condiviso) e il resto del genoma (indicato come strato generale). Abbiamo estratto i p-valori separatamente per ciascuno dei due strati definiti e poi calcolato i q-valori attraverso il quadro FDR. Questi valori q stratificati erano le quantità finali considerate per la significatività statistica (<10-2) per individuare gli SNPs sotto selezione e fonderli nelle regioni genomiche corrispondenti.

Per verificare la presenza di firme convergenti di selezione che differenziano gli animali selvatici da quelli domestici in entrambi i generi, abbiamo esaminato la relazione tra la soglia di significatività applicata ai valori q (che abbiamo fatto variare da 0,2 a 0,002) in un genere e la probabilità stimata che uno SNP sia selezionato nello strato condiviso dell’altro genere utilizzando l’approccio di Storey et al.70 approccio. Un aumento della probabilità dedotta con una diminuzione della soglia applicata al valore q (aumento della severità) indica che più la regione è significativa in un genere, più è probabile trovare SNP significativi nell’altro genere.

Abbiamo filtrato i segnali di selezione che non erano coerenti tra i tre gruppi domestici. Per ogni regione rilevata, abbiamo usato gli aplotipi fasati degli individui che sono stati raggruppati usando alberi Neighbor-Joining basati sulla percentuale di identità tra le sequenze. Solo le regioni che mostravano segnali coerenti sono state mantenute (Fig. 5 supplementare).

Al fine di dedurre se i segnali di selezione rilevati con hapFLK indicavano un rilassamento della selezione o una selezione positiva nei domestici, abbiamo stimato la differenza nella diversità nucleotidica (π) su ogni regione putativa sotto selezione tra i gruppi selvatici e domestici. Abbiamo espresso questa differenza come l’indice Δπ, che è stato calcolato per ogni regione genomica come la differenza tra π calcolato per il gruppo selvatico e la media di π per i gruppi domestici iraniani e marocchini, meno la differenza di π tra questi due gruppi calcolata sull’intero genoma:

Un valore negativo indicherebbe che la diversità nucleotidica è più bassa nel gruppo selvatico rispetto alla media dei due gruppi domestici, e sarebbe considerato come indicativo di un rilassamento della selezione in questi ultimi gruppi, una selezione diversificante nei domestici o una selezione positiva nei selvatici. Al contrario, un valore positivo indicherebbe una selezione positiva direzionale o stabilizzante avvenuta nei gruppi domestici. Abbiamo anche utilizzato il clustering degli aplotipi per verificare manualmente in ogni regione se lo sweep selettivo rilevato confermava le indicazioni date dall’indice Δπ.

Abbiamo condotto interpretazioni funzionali come segue. Per ogni regione sotto selezione, abbiamo considerato la regione più 50 kb su ogni lato per identificare i ruoli funzionali e 5 kb a monte e a valle dei geni e abbiamo valutato la sovrapposizione tra queste coordinate per mantenere i geni di interesse. Infine abbiamo considerato che un gene era legato a una data regione rilevata quando le posizioni della regione e del gene erano sovrapposte. Abbiamo poi valutato quale gene fosse il più probabile bersaglio della selezione considerando il gene più vicino al segnale superiore, cioè la posizione del valore q più basso all’interno della regione. I geni sono stati annotati funzionalmente usando Uniprot (http://www.uniprot.org/), considerando il loro coinvolgimento in 30 termini figli (cioè i discendenti diretti dei termini) della categoria “Processo biologico” (cioè GO:0008150). Abbiamo recuperato tutti i termini GO corrispondenti a ciascun gene (Dati supplementari 4) per 30 delle 33 categorie, perché non abbiamo considerato tre termini che non erano coinvolti nelle funzioni dei mammiferi (cioè, GO:0006791 utilizzazione dello zolfo, GO:0006794 utilizzazione del fosforo, GO:0015976 utilizzazione del carbonio). Abbiamo eseguito due test χ2 per confrontare le distribuzioni dei geni nelle categorie GO, cioè (i) i geni sotto selezione dalle regioni specifiche del genere contro quelli delle regioni omologhe, e (ii) tutti i geni sotto selezione contro i 18.689 geni umani associati ai termini GO in Swiss-Prot. Al fine di interpretare le funzioni dei geni in un contesto zootecnico, abbiamo anche recuperato le informazioni disponibili dalla letteratura sui loro effetti fenotipici.

Infine, per trovare gli SNPs all’interno delle regioni precedentemente individuate che erano i più differenziati tra i gruppi selvatici e domestici, abbiamo usato la statistica FLK. Come per hapFLK, essa rappresenta la deviazione delle frequenze alleliche dei singoli marcatori rispetto al modello neutro stimato dalla matrice di parentela65. La stessa procedura è stata utilizzata per adattare i punteggi delle due analisi a una distribuzione χ2 e combinare i p-valori ottenuti come è stato utilizzato per il test hapFLK. Tuttavia, la distribuzione non uniforme dei p-valori ha precluso l’applicazione del quadro FDR e abbiamo selezionato SNPs all’interno delle regioni rilevate con hapFLK mostrando p-valori <10-4. Per questi SNPs abbiamo usato le annotazioni del Variant Effect Predictor (VEP)71 che sono state generate dall’annotazione del genoma Ensembl v74 delle pecore OARv3.1 per Ovis (http://www.ensembl.org/Ovis_aries/Tools/VEP) e dall’annotazione del genoma della capra CHIR1.0 prodotta dalla pipeline di annotazione del genoma eucariotico NCBI per Capra (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/). Gli SNP sono stati classificati in posizioni intergeniche, a monte e a valle (comprese le UTR), introniche ed esoniche. Le differenze tra le distribuzioni di SNPs con p-valori FLK <10-4 e tutti gli SNPs utilizzati per rilevare le firme di selezione sono state esaminate con un test χ2.

Disponibilità dei dati

Le sequenze e i metadati generati per i 73 campioni Ovis e 72 Capra utilizzati in queste analisi sono disponibili al pubblico. Le informazioni generali e tutti i file vcf possono essere trovati sul sito web di Ensembl (http://projects.ensembl.org/nextgen/). Tutti i file Fastq, i file Bam e gli assemblaggi de novo di O. orientalis e C. aegagrus possono essere trovati sull’European Nucleotide Archive (https://www.ebi.ac.uk/ena) sotto il codice di adesione del progetto Nextgen (PRJEB7436).