Miglioramento del genoma di riferimento di Brassica rapa attraverso il sequenziamento a singola molecola e le tecnologie di cattura della conformazione del cromosoma

Assemblaggio del genoma

Per guidare l’assemblaggio del genoma, abbiamo stimato la dimensione del genoma di B. rapa attraverso la citometria a flusso usando il riso come riferimento. Inizialmente abbiamo stimato che B. rapa ha una dimensione del genoma di 455 Mb (Tabella supplementare S1). Ulteriori indagini che coinvolgono i calcoli per la lunghezza totale della mappa di consenso generata sulla base dei dati BioNano hanno indicato una dimensione del genoma di 442,9 Mb (Tabella supplementare S2). Entrambe le stime erano più piccole della dimensione precedentemente riportata di 52915 o 485 Mb2.

Abbiamo assemblato il genoma di B. rapa utilizzando ~57-fold copertura di subread di sequenziamento PacBio (~25.88 Gb), ~ 456 volte la copertura dei dati BioNano (~ 207,70 Gb), e ~ 164 volte la copertura di Hi-C legge (~ 74,64 Gb). L’assemblea risultante consisteva di 1476 contigs, con un N50 contig di 1,45 Mb e una lunghezza totale di 351,06 Mb (Tabella 1). Successivamente, abbiamo rilevato discrepanze all’interno di 22 contigs utilizzando le letture Hi-C (Tabella supplementare S4). Invece di rimuovere questi contigs, li abbiamo divisi nelle regioni di conflitto; i dati per Contig01464 sono mostrati come esempio (Figura S1 supplementare).

Dopo lo scaffolding e la stima delle dimensioni del gap utilizzando le mappe BioNano e le letture mate-pair (da BRAD, http://brassicadb.org), abbiamo ottenuto 1301 scaffold con uno scaffold N50 di 4,44 Mb (Tabella 1). Per assegnare gli scaffold risultanti alle loro posizioni cromosomiche, abbiamo ancorato questi scaffold usando i dati Hi-C e la mappa genetica migliorata (vedi Metodi). Abbiamo ancorato 298,19 Mb di sequenza su dieci cromosomi che includevano 200 scaffold raggruppati dai dati Hi-C e 8 scaffold assegnati dalla mappa genetica. Il nostro assemblaggio finale, chiamato genoma di B. rapa v3.0, ha totalizzato 353,14 Mb di sequenza con 396 lacune (2,08 Mb) (Tabella 1). Il genoma di B. rapa v3.0 è più lungo di v1.5 ma più corto di v2.5.

Per valutare la qualità del genoma di B. rapa v3.0, abbiamo usato varie fonti di dati. In primo luogo, abbiamo convalidato la completezza del nostro assemblaggio attraverso la ricerca di geni eucarioti di base (CEGs) utilizzando CEGMA16. Un totale di 247 su 248 CEGs erano completi, e 1 CEG era parziale, indicando che tutti i CEGs potrebbe essere rilevato nel nostro assembly (Tabella supplementare S6). Successivamente, la qualità del genoma è stata testata abbinando le sequenze di expressed sequence tags (ESTs) di B. rapa (scaricate da dbEST a NCBI), che ha mostrato che il 99,34% delle ESTs potrebbe essere trovato nel genoma di B. rapa v3.0 appena assemblato.

Miglioramento della contiguità

Il genoma di B. rapa v3.0 ha migliorato la contiguità in termini di gaps e dimensioni contig. Il genoma di B. rapa v1.5 è stato generato da sequenze Illumina, mentre per l’assemblaggio v2.5 sono state utilizzate più letture Illumina e una quantità relativamente piccola di dati di sequenza PacBio. Questi due assemblaggi hanno delle limitazioni dovute alla loro frammentazione e bassa contiguità (Tabella 1). Combinando sequenziamento a singola molecola, mappatura ottica e tecnologia Hi-C, B. rapa genoma v3.0 rappresenta un ~ 27 volte (contig N50: 1446 Kb vs 53 Kb, v2.5) e ~ 31 volte (contig N50: 1446 Kb vs 46 Kb, v1.5) miglioramento della contiguità rispetto ai due precedenti assemblee (Tabella 1). C’erano solo 396 lacune nella v3.0, comprese le lacune di dimensioni note (122 da BioNano e 74 da dati mate-pair) e sconosciute (190 da Hi-C scaffold joining e 10 da genetica map joining). Rispetto agli assemblaggi precedenti, v3.0 ha ~10 volte (5.89 Kb vs. 60.59 Kb, v2.5) e ~7 volte (5.89 Kb vs. 40.09 Kb, v1.5) il miglioramento nella dimensione delle lacune per Mb rispetto ai due assemblaggi precedenti (Tabella 1). In termini di numero di lacune per Mb, v3.0 è superiore a v2.5 e v1.5, rispettivamente, con ~ 23 volte (1.15 vs. 25.98, v2.5) e ~ 35 volte (1.15 vs. 40.09, v1.5) meno lacune per Mb (Tabella 1).

Per valutare la contiguità e la precisione di scaffold ordinamento delle tre versioni della B. rapa, abbiamo prima ricostruito le mappe genetiche basate sui tre assemblaggi utilizzando lo stesso set di dati di risequenziamento di una popolazione aploide raddoppiata (DH) derivata da un incrocio di due linee di cavolo cinese17. Abbiamo poi valutato le posizioni dei binmarkers sulle mappe genetiche integrandole con le corrispondenti mappe fisiche. Dei 892 binmarkers nel nostro assemblaggio, 877 binmarkers (98,3%) sono stati mappati nella mappa genetica. Il nostro assemblaggio concordato con la mappa genetica per 801 binmarkers (91,3%), indicando l’alta qualità della v3.0 (Fig.1; Tabella supplementare S7). Tuttavia, abbiamo notato che 76 (8,7%) binmarkers sui cromosomi A05, A08, e A09 mappato in posizioni ambigue nella mappa genetica. Queste regioni contenevano sequenze ripetute, specialmente nelle regioni centromeriche, come descritto nella seguente analisi. Tuttavia, queste regioni in conflitto sono state coperte dalle letture PacBio e/o dalle mappe BioNano; i dati per il cromosoma A08 nella v3.0 sono mostrati come esempio (Figura supplementare S2).

C’erano 1092 marcatori sulla mappa genetica di v2.5 e 866 marcatori sulla mappa genetica di v1.5. Tuttavia, abbiamo potuto mappare solo l’88,7% dei marcatori (969 di 1092) e il 92,3% dei marcatori (799 di 866) sulla mappa genetica di v2.5 e v1.5, rispettivamente (Tabella supplementare S7). Abbiamo trovato che il 15,1% dei binmarkers (166 di 969) in v2.5 erano discrepanti, compresi 146 binmarkers con distanze genetiche e fisiche disordinate all’interno dello stesso cromosoma (intracromosoma) e 20 binmarkers con distanze genetiche e fisiche incoerenti su diversi cromosomi (intercromosoma) (Figura supplementare S3; Tabella supplementare S7). Per la v1.5, il 10,0% dei binmarkers (80 su 799) erano discrepanti, compresi 71 binmarkers di intrachromosome e 9 binmarkers di interchromosome (Figura supplementare S4; Tabella supplementare S7). Tuttavia, v3.0 conteneva il meno conflittuale intrachromosomal binmarkers (8.7%, 76 di 877 binmarkers) e nessun binmarkers discordante interchromosomal (Tabella supplementare S7), indicando che B. rapa genoma v3.0 ha una maggiore contiguità rispetto alle due assemblee precedenti. Presi insieme, queste convalide indipendenti suggeriscono che B. rapa genoma v3.0 ha la più alta contiguità e il miglior ordinamento di scaffold tra i tre assemblaggi B. rapa.

Confronto di annotazione genoma

Abbiamo predetto e annotato i modelli di gene come precedentemente descritto6. Abbiamo identificato un totale di 45.985 modelli di geni codificanti proteine nella v3.0, che rappresentava il 14,74% dell’assemblaggio del genoma (Tabella 1). Nel nostro assemblaggio, 98.75% (45.411 di 45.985) dei geni sono stati annotati sui cromosomi, e solo 1.25% (574 di 45.985) si trovava su scaffold. I geni annotati de novo nella v3.0 sono stati nominati seguendo lo standard di nomenclatura dei modelli genici per i genomi di riferimento di Brassica (http://www.brassica.info/info/genome_annotation.php). Il numero di modelli genici nel nuovo assemblaggio è superiore a quello della v1.5 (41.020 geni) ma inferiore a quello della v2.5 (48.826 geni) (Tabella 1). Per valutare ulteriormente la qualità dell’annotazione, è stato eseguito un confronto con l’annotazione dei precedenti assemblaggi utilizzando BUSCO18, che si basa su un benchmark di 1440 geni vegetali conservati. Circa il 97,7% di questi geni vegetali conservati sono stati identificati, e 1,7% sono stati rilevati come frammenti presentati in v3.0 (Supplementary TableS11).

Un’analisi di sintonia del genoma è stata eseguita tra i tre assemblee utilizzando SynOrths19 per identificare coppie di geni sintenici e tandem gene array. Un totale di 2077 array tandem (corrispondenti a 4963 geni tandem) sono stati identificati nella v3.0. Lo stesso numero di array tandem (2077 array corrispondenti a 5004 geni) è stato rilevato anche in v1.5. Una valutazione della sinteticità a livello di genoma ha indicato che 1539 matrici in tandem (corrispondenti a 3757 geni) nella v3.0 erano sinteniche a 1494 matrici in tandem (corrispondenti a 3670 geni) nella v1.5. Tuttavia, più array tandem (3535 array, 8002 geni) sono stati identificati nella v2.5 (Tabella 1). Abbiamo rilevato lacune nelle regioni dei geni tandem superflui in v2.5, mentre non sono state trovate lacune né in 3.0 né in v1.5 (Fig. 2a). Queste lacune possono essere il risultato di errori di assemblaggio prodotto da chiusura gap utilizzando PacBio legge in v2.5, che a sua volta ha portato alla annotazione non valida di geni tandem. Per altri geni tandem senza lacune, abbiamo osservato che i singoli geni in v3.0 e v1.5 sono stati annotati come due o più geni in v2.5 (Fig. 2b).

Prendendo ogni array tandem come un singolo locus genico, sono rimasti 43.099 geni nella v3.0, 44.359 geni nella v2.5 e 38.093 geni nella v1.5 (tabella 1). Abbiamo quindi eseguito un’analisi della sinteticità genica, che ha rivelato che 39.858 geni (92,48%) nella v3.0 sono serviti come controparte a 40.442 (91,17%) e 35.464 geni (93,10%) nella v2.5 e v1.5, rispettivamente. Dopo aver confrontato i geni annotati con quelli delle prime versioni, abbiamo identificato 3241 geni specifici della versione nella v3.0 rispetto sia alla v2.5 che alla v1.5. Di questi, 2380 geni sono stati supportati da prove da corrispondenti letture mRNA di B. rapa (da BRAD, http://brassicadb.org/), e 2295 geni sono stati supportati da sequenze di proteine di altre specie di Brassicaceae (Tabella supplementare S12). In totale, 89,10% (2888 di 3214) dei geni specifici della versione in v3.0 sono stati supportati dai dati mRNA di B. rapa o le sequenze proteiche di altre specie di Brassicaceae, mentre solo 10,90% (326 di 3214) dei geni non sono stati supportati.

Un nuovo evento di espansione LTR-RT identificato nel montaggio aggiornato

Abbiamo annotato TEs in v3.0 utilizzando gli stessi metodi come precedentemente riportato20. Un totale di 235.683 TE sono stati identificati da 1244 famiglie nella v3.0, e 562 famiglie TE uniche sono state trovate rispetto alla v2.5 e v1.5. In v3.0, i TE rappresentano 37.51% (134 Mb) del genoma assemblato, che era più alto rispetto ai precedenti assemblaggi (32.30%, 126 Mb, v2.5; 25.44%, 72 Mb, v1.5)2,6. Nel nostro nuovo assemblaggio, i TE più abbondanti sono LTR-RT, che copre una lunghezza totale di 57,64 Mb e rappresenta il 16,32% del genoma assemblato. Non-LTR-RT ripete (LINEs e SINEs) rappresentano 3.10% del nostro assemblaggio (Figura supplementare S5). Abbiamo rilevato trasposoni di DNA corrispondenti a 26,35 Mb, che costituiscono il 7,46% del genoma assemblato (Figura supplementare S5). Un elenco completo di TEs identificati e ripetizioni in v3.0 può essere trovato nella tabella supplementare S13. Inoltre, abbiamo identificato un totale di 1231 miRNA, 1281 tRNA, 2865 rRNA e 3737 snRNA nel genoma di B. rapa v3.0 (Tabella supplementare S19).

Nel nostro attuale assemblaggio, abbiamo annotato più LTR-RT (57 Mb) rispetto alla v2.5 (44 Mb) e v1.5 (18 Mb). Abbiamo identificato 51.062 LTR-RT non intatte nella v3.0. Ulteriori analisi hanno rivelato che il 65,27% (33.672 di 51.602) delle LTR non intatte si trovavano sui dieci cromosomi, mentre il 34,73% (17.922 di 51.602) delle LTR-RT non intatte sono state trovate sulle impalcature non ancorate. Usando lo stesso metodo6, un totale di 13.318 LTR-RT intatte sono state annotate nella v3.0. Tuttavia, c’erano solo 4129 e 801 intact LTR-RTs in v2.5 e v1.5, rispettivamente6. Ulteriori analisi hanno rivelato che solo il 18,19% delle LTR-RT intatte (2423 di 13.318) si trovavano sui dieci cromosomi, mentre la maggior parte (81,81%, 10.895 di 13.318) LTR-RT intatte sono state trovate sulle impalcature non ancorate in v3.0.Il tempo di inserimento delle LTR-RT intatte è stato calcolato come precedentemente descritto4, che ha indicato che il genoma di B. rapa ha subito tre ondate di espansione LTR-RT da quando si è differenziato da B. oleracea (Fig. 3). Queste LTR-RT intatte avevano un’età media di inserimento di 1,88 milioni di anni fa (MYA), con un’età mediana di inserimento di 1,59 MYA. Inoltre, abbiamo trovato più intatto LTR-RTs con lunghezze diverse in v3.0 rispetto a in v2.5 e v1.5 (Figura supplementare S6).

Con questi intatti LTR-RTs, un nuovo evento di espansione LTR-RT è stato identificato nel genoma B. rapa. Abbiamo designato 3155 eventi di inserzione LTR-RT intatti da 0 MYA a 0,4 MYA come “giovane espansione” con una lunghezza media di 8135 bp e una data media di inserimento di 0,2 MYA; 2283 eventi di inserzione LTR-RT intatti da 1,0 MYA a 1,4 MYA come “espansione media” con una lunghezza media di 11.902 bp e una data media di inserimento di 1,2 MYA; e 1444 eventi di inserzione LTR-RT intatti da 3.0 MYA a 3,4 MYA come “espansione antica” con una lunghezza media di 9823 bp e data di inserimento (Fig. 3).Le espansioni giovani e antiche corrispondono strettamente agli eventi di espansione precedentemente identificati; l’espansione media è stata identificata per la prima volta nel genoma di B. rapa e ha un tempo di inserimento simile a quello dell’evento di espansione LTR-RT intatto in B. oleracea. Inoltre, 1778 Ty1/Copia-like LTR-RTs e 4179 Ty3/Gypsy-like LTR-RTs sono stati identificati in v3.0, che è molto più di quelli identificati negli assemblaggi precedenti (353 Ty1/Copia e 632 Ty3/Gypsy in v2.5, 260 Ty1/Copia e 162 Ty3/Gypsy in v1.5) (Tabella supplementare S20; Figura supplementare S7, S8). In generale, c’erano più Ty3/Gypsy-like LTR-RTs che Ty1/Copia-like LTR-RTs (Tabella supplementare S20). Rispetto a v2.5 e v1.5, Ty3/Gypsy-like LTR-RTs in v3.0 erano ovviamente aumentati dal 5 MYA (Figura supplementare S7), mentre Ty1/Copia-like LTR-RTs erano aumentati dal 2.2 MYA (Figura supplementare S8). Dagli alberi filogenetici, abbiamo trovato che ogni gruppo di LTR-RTs aveva più copie in v3.0 che in v2.5 e v1.5 (Tabella supplementare S21, S22; Figura supplementare S9, S10).

Bocchi di genoma e centromeri nel genoma di B. rapa

Abbiamo studiato le relazioni dei blocchi di genoma usando l’assemblaggio aggiornato v3.0. Per definire i blocchi del genoma e i centromeri nel genoma di B. rapa v3.0, abbiamo prima costruito i tre sottogenomi (LF, MF1 e MF2) sulla base della relazione sintenica tra v3.0 e A. thaliana (Figura supplementare S11; Tabella supplementare S14). Abbiamo rilevato 71 dei 72 (3 × 24) blocchi genomici previsti in v3.0, e la maggior parte di essi sono stati disposti in linea con quelli precedentemente riportati in rif.10 (Fig. 4; Tabella supplementare S15). In v3.0, i due nuovi blocchi genomici frammentati F (LF) e F (MF1) sono stati identificati sui cromosomi A01 e A05 e non sono stati osservati in rif.10. Non abbiamo potuto rilevare due blocchi di genoma molto piccoli precedentemente descritti nella v3.0, il blocco C (MF2) sul cromosoma A07 e il blocco B (MF1) sul cromosoma A08 (rif.10). Tuttavia, nel nostro assemblaggio, i blocchi di genoma N/M (MF1), O/P (LF) e A/C(LF) sono stati disposti rispettivamente sui cromosomi A01, A09 e A10, mentre erano ordinati su lati opposti in rif.10. I tre piccoli blocchi di genoma adiacenti (S (MF2), T (MF2) e B (MF1)) sul cromosoma A08 della v3.0 erano ordinati S/T/B, mentre questi erano disposti come T/B/S in rif.10.

Abbiamo anche confrontato l’orientamento dei blocchi del genoma in v3.0 con quello in rif.10. I blocchi del genoma W (MF2) sul cromosoma A02, così come G (LF) ed E (LF) sul cromosoma A07, sono risultati invertiti rispetto agli altri blocchi che hanno avuto origine da un singolo cromosoma ACK. Tuttavia, l’orientamento del blocco del genoma P (LF) sul cromosoma A09 e tre blocchi di V in v3.0 erano in direzione avanti, mentre questi erano invertiti in rif.10. Questi risultati sono stati ulteriormente supportati dalle mappe genetiche di v3.0 e v1.5, rispettivamente.

Abbiamo determinato con precisione la posizione dei centromeri di tutti i cromosomi in v3.0. Selezionando le sequenze di ripetizioni centromeriche precedentemente determinate, comprese le ripetizioni centromeriche satellite CentBr, CRB, TR238, e PCRBr21,22,23, abbiamo identificato i segnali per tutte le 21 regioni paleocentromeriche in v3.0, mentre tre regioni paleocentromeriche non sono state rilevate in rif.10 (Fig. 4, Tabella supplementare S16). L’analisi dei paleocentromeri ha indicato che i dieci centromeri estanti di B. rapa sono stati tutti ereditati dai 21 paleocentromeri. In v3.0, i centromeri dei cromosomi A01, A03, A04, A05, A06, A07 e A10 avevano gli stessi blocchi di genoma associati che fiancheggiano i centromeri corrispondenti come riportato in rif.10 (Fig. 4). Tuttavia, il centromero sul cromosoma A02 si trovava tra i blocchi di genoma P (MF2) e V (MF1), e il centromero sul cromosoma A09 era situato tra i blocchi di genoma P (LF) e B (LF), mentre questi erano considerati paleocentromeri in rif.10 (Fig. 4). Il centromero sul cromosoma A08 si trovava tra i blocchi di genoma T (MF2) e B (MF1), piuttosto che tra i blocchi di genoma C (MF1) e T (MF2), come riportato in rif.10. Inoltre, c’erano 1188 geni rilevati all’interno di regioni centromeriche in v3.0, mentre solo 740 geni sono stati rilevati in rif.10 (Tabella supplementare S17).

Per valutare il nostro assemblaggio per quanto riguarda i centromeri in v3.0, abbiamo analizzato le caratteristiche di sequenza delle regioni centromeriche. Abbiamo trovato che un numero significativamente più alto di TE e di ripetizioni specifiche del centromero sono state mappate nelle regioni centromeriche rispetto ad altre parti dei cromosomi, e la densità genica e il tasso di ricombinazione erano nettamente inferiori alle regioni centromeriche annotate nella v3.0 (Fig. 5). Inoltre, sono state rilevate più ripetizioni centromeriche nelle regioni centromeriche nella v3.0 rispetto a quelle riportate nel rif. 10 (Tabella supplementare S17).