Heading

Stima e partizionamento basato sull’età dei livelli genomici individuali di inbreeding nel bestiame blu belga

Siamo partiti utilizzando un modello Mix14R (con R k che varia da 2 a 819).di partizione basata sull’età dei livelli individuali di inbreeding genomico nella popolazione di bovini Belgian Blue Beef

Abbiamo iniziato utilizzando un modello Mix14R (con R k che va da 2 a 8192) per stimare la proporzione del genoma appartenente a diverse classi HBD per i 634 tori BBB (Fig. 1a), che permette la stima del coefficiente di inbreeding rispetto alle diverse popolazioni di base come spiegato nei Metodi (Fig. 1b). Considerando tutte le classi HBD, la frazione del genoma che è HBD (corrispondente al coefficiente di inbreeding stimato con la popolazione di base più lontana) è stata pari a 0,303 in media (variando da 0,258 a 0,375), con un importante contributo delle classi HBD con alti tassi R k (R k > 256) che rappresentano in media il 71,4% della proporzione totale HBD. Questi piccoli ROH riflettono la storia della popolazione (inbreeding di fondo e linkage disequilibrium associati alla dimensione effettiva della popolazione (Ne) del passato) meglio della variazione individuale. Le classi associate a tassi R k più piccoli (cioè con segmenti HBD più lunghi) hanno rappresentato una proporzione minore della proporzione HBD totale (il coefficiente medio di inbreeding era pari a 0,054 e 0,087 quando si includevano le classi HBD con R k ≤ 32 e R k ≤ 256, rispettivamente, e si impostava la popolazione di base circa 16 o 128 generazioni fa) ma presentavano più variazione tra gli individui. Per esempio, il coefficiente di inbreeding associato ad antenati comuni risalenti fino a circa quattro generazioni fa (corrispondenti a classi HBD con R k ≤ 8) variava da 0,000 a 0,137. Per i tori nati dal 1980 al 2010, la percentuale del genoma in segmenti HBD è aumentata del 3,3% (+ 0,11% per anno), cioè approssimativamente dal 28 al 31% (vedi file aggiuntivo 1: Fig. S1a). Tuttavia, la tendenza per le classi HBD più recenti (R k ≤ 32) è stata più pronunciata (vedi file aggiuntivo 1: Fig. S1b), cioè da quasi 0 al 6% (+ 0,20% all’anno) e corrispondeva più strettamente alla tendenza osservata con i coefficienti di consanguineità basati sul pedigree (vedi file aggiuntivo 1: Fig. S1c). I tori nati prima del 1980 presentavano poche prove di autozigosi recente rispetto ai tori moderni.

Per valutare il contributo di ogni classe HBD alla percentuale del genoma in segmenti HBD e alla sua variazione in BBB cattle, abbiamo diviso la frazione totale del genoma in classi HBD in quattro classi principali (classi HBD molto recenti con R k = da 2 a 8, classi HBD recenti con R k = da 16 a 64, classi HBD antiche con R k = da 128 a 512, e classi HBD molto antiche con R k = da 1024 a 8192), con ogni gruppo avente tre classi HBD tranne l’ultimo con quattro classi HBD. La frazione media del genoma associata a ciascuna di queste classi principali (ordinate da recenti ad antiche) era pari a 0,027 (SD = 0,029), 0,041 (SD = 0,019), 0,054 (SD = 0,013) e 0,180 (SD = 0,011). Si noti che alte proporzioni di segmenti HBD molto recenti sono meccanicamente associate a proporzioni inferiori di segmenti HBD molto antichi (r = – 0,407) perché i segmenti HBD recenti mascherano segmenti HBD più antichi. Sebbene la percentuale del genoma in classi HBD associate ad antenati comuni recenti rappresenti solo il 22,6% dell’autozigosità totale, essa mostra più variazione individuale di quella in classi più antiche (più del 50% della varianza totale è associata a classi HBD molto recenti). Adattando un modello lineare, abbiamo stimato che le classi HBD molto recenti rappresentano il 59% della variazione totale dell’autozigosità e che l’aggiunta di classi HBD recenti al modello aumenta questo valore all’83%. Allo stesso modo, le correlazioni tra i coefficienti di inbreeding misurati rispetto a diverse popolazioni di base (ad esempio, includendo diverse classi HBD nel calcolo) con i coefficienti di inbreeding stimati utilizzando tutte le classi HBD sono aumentati bruscamente da 0,16 per le stime basate sulla prima classe (R k = 2) a 0,77 per i coefficienti di inbreeding stimati includendo le classi HBD con un R k ≤ 8 e a 0,90 con un R k ≤ 16, per poi migliorare solo marginalmente aggiungendo più classi HBD (Fig. 2). La diminuzione della correlazione osservata a R k = 1024 deriva dal fatto che l’autozigosi antica si concentra a R k = 1024 per alcuni individui e a R k = 2048 per altri.

Confronto dei risultati per i bovini BBB con quelli di altre razze

Per determinare se livelli e modelli comparabili di autozigosi sono osservati anche in altre razze di origine europea, abbiamo applicato lo stesso modello a 10 razze genotipizzate con lo stesso array (Fig. 3). Nella maggior parte di queste razze, i coefficienti di inbreeding stimati rispetto alle diverse popolazioni di base sono aumentati moderatamente fino a FG-256 (ad esempio, classe HBD con R k ≤ 256 inclusa nella stima) e più fortemente con le popolazioni di base più vecchie (da FG-512 a FG-2048), che includono molte più generazioni di antenati. Grandi differenze nei coefficienti di inbreeding sono state osservate con popolazioni di base relativamente recenti (FG-64, circa 32 generazioni fa), che vanno da 0,013 e 0,042 in Piemontese e Limousin a 0,164 e 0,200 in Jersey e Hereford. Alcuni individui Hereford hanno presentato coefficienti di inbreeding estremi stimati con popolazioni di base recenti (vedi file aggiuntivo 2), cioè fino al 40% per FG-8 (ad esempio, circa quattro generazioni indietro). Parte degli individui Hereford di questo set di dati provengono dalla linea Hereford 1, una linea inbred, il che indica che il nostro modello cattura correttamente gli eventi estremi ma anche che gli individui genotipizzati inclusi in questo studio non sono necessariamente rappresentativi della razza.

Stima dei coefficienti di consanguineità e delle probabilità di HBD con diverse densità di SNP

Abbiamo applicato un modello Mix14R utilizzando diverse densità di SNP, cioè.e., da LD (6844 SNPs) a HD (601,226 SNPs) sul dataset 634 BBB e anche a WGS (5,653,911 SNPs) per i 50 individui sequenziati dell’intero genoma. I coefficienti medi di inbreeding stimati misurati rispetto alle diverse popolazioni di base (Fig. 4) e Additional file 3: Fig. S2 erano simili tra i pannelli SNP per le popolazioni di base più recenti (FG-32). Per le popolazioni di base più antiche, meno autozigosi sono state catturate con il pannello LD con differenze marcate per le classi HBD antiche che sono state catturate solo con i pannelli HD o WGS. Una tendenza simile è stata osservata con il pannello 50 K, ma i coefficienti di inbreeding medi erano simili a quelli del pannello HD fino a FG-256 (circa 128 generazioni indietro). I coefficienti di inbreeding medi stimati utilizzando la popolazione di base più remota e i pannelli LD, 50 K e HD erano pari a 0,060, 0,093 e 0,303, rispettivamente (quando stimati sui soli 50 individui sequenziati, questi valori erano pari a 0,047, 0,101 e 0,309, rispettivamente, e a 0,359 con il pannello WGS). La popolazione di base è quindi una funzione dei più piccoli segmenti HBD che possono essere catturati dal pannello utilizzato. Le correlazioni tra questi coefficienti di inbreeding stimati con i diversi panel sono state elevate, cioè 0,934 (LD-HD), 0,944 (LD-50 K) e 0,975 (50 K-HD). Nonostante i coefficienti di inbreeding molto più bassi ottenuti con il pannello 50 K, esso cattura essenzialmente tutta la variazione individuale ottenuta con un pannello HD, in accordo con la precedente osservazione che la maggior parte della variazione era associata alle classi HBD recenti.

Abbiamo quindi utilizzato l’algoritmo Viterbi per identificare i segmenti HBD con diversi pannelli SNP (Tabella 1). L’algoritmo di Viterbi classifica ogni posizione SNP come HBD o non-HBD mentre l’algoritmo forward-backward fornisce la probabilità locale di HBD. Come ci si aspettava, un maggior numero di segmenti HBD e più corti sono catturati con pannelli a più alta densità. Con il pannello HD, è stata catturata una proporzione limitata di segmenti estremamente piccoli (pochi kb). La lunghezza della maggior parte dei segmenti variava da 10 a 500 kb, con più della metà più corta di 100 kb, ma tali segmenti non hanno necessariamente il maggior contributo alla percentuale totale del genoma in classi HBD poiché le classi con meno segmenti ma più lunghi possono rappresentare una grande proporzione di autozigosi. Abbiamo anche osservato segmenti HBD estremamente lunghi (> 50 Mb), che hanno confermato la presenza di autozigosi recenti (il segmento HBD più lungo era lungo più di 90 Mb). In media, ognuno dei 634 tori aveva 4,25 segmenti HBD più lunghi di 10 Mb e associati a un antenato comune presente circa cinque generazioni prima. Il numero di tali segmenti HBD variava da 0 a 14 per individuo. Sessantuno tori avevano anche uno o più (fino a tre) segmenti HBD più lunghi di 50 Mb. Con i pannelli 50 K e LD, più del 99% dei segmenti identificati erano più lunghi di 100 e 500 kb, rispettivamente (con un picco nelle classi da 1 a 5 Mb e da 5 a 10 Mb, rispettivamente), e solo una frazione dei segmenti è stata catturata rispetto a quando è stato utilizzato il pannello HD. In particolare, la maggior parte dei segmenti HBD più corti di 1 Mb non sono stati identificati. A basse densità di SNP, i segmenti più piccoli semplicemente non vengono catturati perché non contengono alcun SNP o sono troppo pochi. I segmenti di dimensioni intermedie potrebbero non raggiungere alte probabilità di HBD a causa di un minor numero di SNPs nel segmento. Al contrario, la lunghezza di alcuni segmenti HBD può essere sovrastimata quando si usa il pannello LD, per esempio quando non ci sono abbastanza SNPs per identificare piccoli segmenti non-HBD che affiancano segmenti HBD. La figura 5a illustra l’identificazione dei segmenti HBD per un cromosoma. Mostra che (1) più segmenti sono stati identificati ad alta densità, (2) le probabilità di HBD erano più alte con mappe più dense, (3) l’algoritmo Viterbi ha dichiarato alcune posizioni SNP come HBD anche se avevano solo moderate probabilità di HBD, e (4) i confini dei segmenti HBD variavano con la densità del pannello. Allo stesso modo, la Fig. 5b rappresenta i segmenti HBD che sono stati identificati sul cromosoma Bos taurus (BTA) 5 per 50 individui con l’algoritmo di Viterbi con diverse densità di SNP. I risultati sono in accordo con quelli riportati nella tabella 1. Proporzioni maggiori del genoma sono state dichiarate HBD con il pannello HD e i piccoli segmenti HBD hanno rappresentato la maggior parte della differenza con i risultati dei pannelli a densità inferiore. Tuttavia, abbiamo osservato che alcuni segmenti HBD di pochi Mb di lunghezza non sono stati identificati con densità SNP inferiori (e ancora di più con il pannello LD). Come per la Fig. 5a, la lunghezza di alcuni segmenti HBD è sovrastimata quando è stato usato il pannello LD. Abbiamo anche confrontato le probabilità locali di HBD stimate utilizzando il pannello LD o il pannello 50 K con le classi locali di HBD dedotte utilizzando il pannello HD e l’algoritmo Viterbi (Fig. 6). Le probabilità di HBD erano alte per le classi HBD recenti e scendevano per gli antenati comuni più remoti. Come previsto, il pannello LD era efficiente solo per gli antenati comuni più recenti (la probabilità di HBD era 0,90 o superiore quando R k < 16 e ~ 0,50 per R k = 32) mentre il pannello 50 K ha permesso la cattura di autozigosi più antiche (la probabilità di HBD era 0,90 o superiore quando R k < 64 e ~ 0,50 per R k = 128). Altri risultati riguardanti l’età (o la lunghezza) dei segmenti HBD che possono essere catturati con diverse densità di SNP sono descritti in Druet e Gautier .

Confronto dei modelli

Modelli che stimano i tassi R k delle classi HBD (modelli KR)

Per le diverse densità SNP testate e per ogni individuo, abbiamo usato il BIC (vedi ) per selezionare il modello KR con il miglior supporto statistico (cioè.e., con il numero ottimale di classi K, con K – 1 classi HBD e una classe non-HBD) dopo aver stimato il tasso(i) R k per ogni individuo con ogni modello testato. Per ogni pannello SNP, la tabella 2 mostra il numero di volte che un modello è stato selezionato come il migliore per l’individuo analizzato. All’aumentare della densità degli SNP, più generazioni passate possono essere esplorate e il K ottimale aumenta di conseguenza. Nella maggior parte dei casi, i modelli con una classe HBD sono preferiti per il pannello LD, i modelli con due classi HBD per il pannello 50 K, i modelli con tre classi HBD per i pannelli HD e WGS (sebbene il modello con quattro classi HBD sia spesso selezionato anche per quest’ultimo, cioè per 23 individui su 50). Con questi modelli ottimali, la prima classe HBD cattura l’autozigosi più recente (R k da 15 a 20), la seconda classe HBD cattura l’autozigosi che è associata ad antenati comuni da alcune centinaia di generazioni indietro e le classi successive sono associate a R k più elevati (> 1000) (Tabella 2). Le correlazioni dei coefficienti di inbreeding stimati con questi modelli KR selezionati con quelli ottenuti con il modello completo Mix14R (che vanno da 0,981 a 1,000) e il confronto dei coefficienti di inbreeding medi stimati indicano che con questi modelli KR ridotti, possiamo catturare efficacemente l’autozigosità su tutto il genoma. Con modelli 1R e densità SNP basse o moderate, abbiamo osservato una leggera sottostima dei coefficienti di inbreeding rispetto al modello Mix14R e correlazioni leggermente inferiori (ancora sopra 0,98). I tassi R k stimati per ogni individuo con questi pannelli hanno un valore mediano inferiore (rispettivamente 15 e 41 con i pannelli LD e 50 K) rispetto ai tassi R k stimati con i pannelli a densità maggiore (mediana R k > 1000) per i quali il contributo dei ROH più piccoli è molto maggiore. Di conseguenza, alcuni piccoli frammenti non sono stati catturati dal modello a bassa densità, mentre a densità maggiore, i coefficienti di inbreeding sono quasi identici alle stime ottenute con il modello Mix14R. I modelli contenenti due o più classi HBD hanno catturato la stessa quantità di autozigosi del modello Mix14R, indipendentemente dalla densità degli SNP. Anche se il coefficiente di inbreeding è correttamente stimato con un modello 1R (una classe HBD e una non HBD con lo stesso tasso) con i dati WGS, i segmenti HBD identificati tendono ad essere più piccoli poiché i tassi R k stimati sono più alti (cioè, minori lunghezze attese dei frammenti) come mostrato nel file aggiuntivo 4: Fig. S3. Infatti, il modello 1R produce più segmenti lunghi da 10 a 100 kb rispetto al modello Mix14R, ma meno segmenti più lunghi di 100 kb. Quindi, con un modello 1R, i lunghi segmenti HBD potrebbero essere tagliati in frammenti più piccoli in presenza di SNP eterozigoti (forse errori di sequenziamento), mentre con i modelli che includono le classi HBD associate a recenti antenati comuni (con piccoli tassi R k), questi segmenti HBD sono identificati come un unico frammento lungo e recente (perché la pena per terminare e iniziare un nuovo segmento è più alta). La figura 7 illustra questo con un esempio. Infatti, abbiamo osservato un lungo segmento con alte probabilità di HBD anche se ci sono più posizioni in cui la probabilità del genotipo eterozigote non è nulla (ma questo è limitato rispetto alle regioni laterali). Con il modello Mix14R, questo viene considerato come un lungo segmento e la probabilità locale di HBD rimane superiore a 0,99 per l’intera regione (ad eccezione di una regione con cinque SNP eterozigoti consecutivi). Con il modello 1R, le probabilità di HBD scendono ripetutamente a causa di questi SNP eterozigoti e il segmento HBD più lungo viene tagliato in diversi frammenti più piccoli (in base ai risultati dell’algoritmo di Viterbi). Si noti che con il pannello HD, questo individuo è omozigote per tutti i 13.009 SNPs che sono inclusi in questo segmento lungo 56,1 Mb. Come in Fig. 5, notiamo che l’algoritmo di Viterbi classifica alcune posizioni con una bassa probabilità stimata di HBD come HBD.

Modelli con tassi R k predefiniti di classi HBD (modelli MixKR)

Rispetto ai modelli KR, I modelli MixKR presentano il vantaggio di utilizzare le stesse classi HBD per tutti gli individui (i tassi R k delle classi HBD non sono stimati individualmente ma predefiniti dall’utente) e rendono i confronti tra individui più facili (per esempio, confrontare due individui con una singola classe HBD ma con R k = 8 per il primo e R k = 96 per il secondo non sarebbe facile – gli R k stimati vanno da 4 a 1000). Diversi di questi modelli MixKR (con K = 2, 3 e 4) sono stati testati con il pannello LD (Tabella 3) per valutare se i modelli ridotti con tassi predefiniti di classi HBD sono efficienti. Per selezionare questi tassi predefiniti, o abbiamo usato le mediane dei tassi stimati ottenuti da modelli con lo stesso numero di classi (vedi sezione precedente) o abbiamo selezionato alcune classi dal modello MixKR in modo da coprire la gamma di valori stimati (ad esempio, una classe per i segmenti HBD recenti e una per i segmenti HBD antichi). In accordo con le precedenti osservazioni sui modelli KR, i confronti dei coefficienti di inbreeding stimati con quelli ottenuti con il modello Mix14R indicano che i modelli con una sola classe HBD sottostimano leggermente i coefficienti di inbreeding e risultano in correlazioni più basse (> 0,96) rispetto ai modelli con due o più classi HBD (> 0,99). La presenza di più classi HBD (> 2) permette una migliore valutazione dei contributi delle diverse generazioni passate (es, R k = 8 vs 64) ma non fornisce stime migliori del coefficiente di inbreeding a livello genomico.

Confronto con altri stimatori di coefficienti di inbreeding

Mezzi e range dei coefficienti di inbreeding stimati con diversi metodi e con il pannello HD sono in Tabella 4 e le loro correlazioni sono in Tabella 5, e nel file aggiuntivo 5: Tabelle S1 e S2 per gli altri pannelli. Similmente al nostro modello, i modelli basati sull’omozigosi osservata e sul ROH hanno portato ad alti coefficienti di inbreeding (rispettivamente, 0,644 e 0,151 in media) mentre altri stimatori genomici hanno portato a coefficienti di inbreeding centrati intorno a 0 e comprendenti valori negativi. Va notato che valori più alti si ottengono in media (0,268) quando si usano regole meno stringenti per identificare i ROH (per esempio, finestre di 20 SNPs e almeno 10 SNPs per ROH). Abbiamo osservato correlazioni molto elevate tra le stime basate su HMM ed entrambe le misure basate sull’omozigosi (r = 0,95 con FHOM e FExHOM, queste due misure presentano una correlazione di 1 e sono essenzialmente uguali) o su ROH (r = 0,95 con FROH), il che suggerisce che con un gran numero di SNPs, le euristiche semplici (ignorando le frequenze alleliche, la spaziatura degli SNP, ecc.) sono efficienti (FHOM e FROH sono altamente correlati, r = 0,97). La correlazione tra FHOM stimato con i pannelli LD e HD è pari a 0,890, che è leggermente inferiore alla correlazione tra le stime ottenute con l’HMM per questi due pannelli (r = 0,934), il che indica che gli stimatori globali funzionano ancora correttamente con 6844 SNPs in questa popolazione. I metodi ROH basati su regole sono meno efficienti a basse densità di SNP poiché catturano solo i frammenti più lunghi (5 Mb o più e 20 Mb in media) con i parametri utilizzati nello studio attuale (la dimensione predefinita delle finestre in Plink). Infatti, gli stimatori basati su ROH sono raramente utilizzati con il pannello LD nel bestiame, anche se più segmenti HBD potrebbero essere identificati con regole meno rigorose, a scapito di un maggior tasso di falsi positivi. A bassa densità di SNP, il quadro HMM fornisce ancora probabilità HBD globali e locali corrette, anche se i segmenti HBD non vengono identificati senza ambiguità.

Le correlazioni delle stime del GRM tradizionale con le nostre stime sono moderatamente alte (r = 0.73) e più basse con stimatori di omozigosi (r = 0,63) e stimatori basati su ROH (0,61). La FGRM è stata calcolata con la formula proposta da , che divide tutti i contributi SNP per lo stesso peso. Quando è stato stimato con la formula alternativa, che divide ogni contributo SNP per il proprio peso 2f i (1 – f i ) (f i è la frequenza di SNP i) come in Amin et al. , le correlazioni erano più basse (cioè, 0,48 con FG, 0,34 con FHOM e 0,33 con FROH). Lo stimatore basato sulle correlazioni unificate tra gameti proposto da Yang et al. ha presentato correlazioni relativamente alte sia con FG che con FGRM (rispettivamente, 0,90 e 0,92) e correlazioni leggermente inferiori con gli altri stimatori (r = 0,87 e 0,85 con FHOM e FROH, rispettivamente).

Le correlazioni di queste stime con i coefficienti di inbreeding del pedigree (considerando solo gli individui nati dopo il 1999 per aumentare la profondità del pedigree) sono anche nella tabella 5. Le correlazioni complessive sono state moderate, con i valori più alti per le correlazioni con le misure basate su omozigosi e ROH (0,55 per entrambe le misure) e valori leggermente inferiori per quelle con lo stimatore basato su HMM (0,46), mentre abbiamo osservato una bassa relazione con FGRM (0,29) e una moderata correlazione con FUNI (0,45). Abbiamo anche confrontato il FPED e i coefficienti di inbreeding stimati con il nostro modello rispetto a diverse popolazioni di base (Fig. 8) e abbiamo scoperto che le correlazioni sono aumentate fino a FG-32 (catturando l’inbreeding dagli antenati circa 16 generazioni indietro) e poi hanno presentato un plateau da FG-32 a FG-256 raggiungendo un massimo a r = 0,56 (cioè, leggermente meglio degli stimatori basati sull’omozigosi). Questa tendenza era attesa poiché FPED è stimato per un numero limitato di generazioni indietro nel tempo. Il numero medio equivalente di generazioni conosciute stimato con PEDIG è stato di 6,3 per i tori nati dopo il 1999 (è aumentato da 5,5 per i tori nati nel 2000 a 7,5 per quelli nati nel 2010) corrispondente in media a FG-16. L’aggiunta della classe HBD R k = 32 permette di catturare i contributi di alcuni rami più vecchi del pedigree e i più piccoli segmenti HBD ereditati da antenati comuni nel pedigree.