Improved Brassica rapa reference genome by single-molecule sequencing and chromosome conformation capture technologies

Genome assembly

Aby pokierować montażem genomu, oszacowaliśmy rozmiar genomu B. rapa za pomocą cytometrii przepływowej, używając ryżu jako odniesienia. Wstępnie oszacowaliśmy, że B. rapa ma genom o wielkości 455 Mb (Tabela Uzupełniająca S1). Dalsze badania obejmujące obliczenia dla całkowitej długości mapy konsensusu wygenerowanej na podstawie danych BioNano wskazały na rozmiar genomu 442,9 Mb (Tabela S2). Oba oszacowania były mniejsze niż wcześniej zgłoszony rozmiar 52915 lub 485 Mb2.

Złożyliśmy genom B. rapa przy użyciu ~57-krotnego pokrycia subreadów sekwencjonowania PacBio (~25.88 Gb), ~456-krotnego pokrycia danych BioNano (~207.70 Gb) i ~164-krotnego pokrycia odczytów Hi-C (~74.64 Gb).Wynikowe złożenie składało się z 1476 kontigów, z N50 kontigu 1.45 Mb i całkowitej długości 351.06 Mb (Tabela 1). Następnie wykryliśmy rozbieżności w obrębie 22 kontigów przy użyciu odczytów Hi-C (Supplementary Table S4). Zamiast usuwać te kontigi, podzieliliśmy je w regionach konfliktu; dane dla Contig01464 są pokazane jako przykład (Supplementary Figure S1).

Po wykonaniu rusztowań i oszacowaniu wielkości luk przy użyciu map BioNano i odczytów z par mate (z BRAD, http://brassicadb.org), uzyskaliśmy 1301 rusztowań o scaffold N50 równym 4,44 Mb (Tabela 1). Aby przypisać powstałe rusztowania do ich pozycji chromosomalnych, zakotwiczyliśmy je przy użyciu danych Hi-C i ulepszonej mapy genetycznej (patrz Metody). Zakotwiczyliśmy 298,19 Mb sekwencji na dziesięciu chromosomach, w tym 200 rusztowań zgrupowanych na podstawie danych Hi-C i 8 rusztowań przypisanych na podstawie mapy genetycznej. Nasza ostateczna asymilacja, nazwana genomem B. rapa v3.0, liczyła 353,14 Mb sekwencji z 396 lukami (2,08 Mb) (Tabela 1). Genom B. rapa v3.0 jest dłuższy niż v1.5, ale krótszy niż v2.5.

Aby ocenić jakość genomu B. rapa v3.0, wykorzystaliśmy różne źródła danych. Po pierwsze, sprawdziliśmy kompletność naszego złożenia poprzez wyszukiwanie rdzeniowych genów eukariotycznych (CEGs) przy użyciu CEGMA16. W sumie 247 z 248 CEG było kompletnych, a 1 CEG był częściowy, co wskazuje, że wszystkie CEG mogły być wykryte w naszej asemblacji (Supplementary Table S6). Następnie sprawdzono jakość genomu poprzez dopasowanie sekwencji znaczników sekwencji wyrażonych (EST) B. rapa (pobranych z dbEST w NCBI), co wykazało, że 99,34% EST można znaleźć w nowo zmontowanym genomie B. rapa v3.0.

Poprawa przyległości

Genom B. rapa v3.0 poprawił przyległość pod względem luk i rozmiarów kontigów. Genom B. rapa v1.5 został wygenerowany z sekwencji Illumina, podczas gdy do złożenia v2.5 użyto więcej odczytów Illumina i stosunkowo niewielką ilość danych sekwencyjnych PacBio. Te dwie asemblacje mają ograniczenia wynikające z ich fragmentacji i niskiej przyległości (Tabela 1). Dzięki połączeniu sekwencjonowania jednocząsteczkowego, mapowania optycznego i technologii Hi-C, genom B. rapa v3.0 reprezentuje ~27-krotną (kontig N50: 1446 Kb vs. 53 Kb, v2.5) i ~31-krotną (kontig N50: 1446 Kb vs. 46 Kb, v1.5) poprawę przyległości w stosunku do dwóch poprzednich złożeń (Tabela 1). Było tylko 396 luk w v3.0, w tym luki o znanych (122 z BioNano i 74 z danych mate-pair) i nieznanych rozmiarach (190 z łączenia rusztowań Hi-C i 10 z łączenia map genetycznych). W porównaniu z poprzednimi asemblacjami, v3.0 ma ~10-krotną (5,89 Kb vs. 60,59 Kb, v2.5) i ~7-krotną (5,89 Kb vs. 40,09 Kb, v1.5) poprawę w wielkości luk na Mb w stosunku do dwóch poprzednich asemblacji (Tabela 1). Pod względem liczby luk na Mb, v3.0 przewyższa odpowiednio v2.5 i v1.5, z ~23-krotnie (1.15 vs. 25.98, v2.5) i ~35-krotnie (1.15 vs. 40.09, v1.5) mniejszą liczbą luk na Mb (Tabela 1).

Aby ocenić przyległość i dokładność uporządkowania rusztowań trzech wersji genomu referencyjnego B. rapa, najpierw zrekonstruowaliśmy mapy genetyczne w oparciu o te trzy złożenia, wykorzystując ten sam zestaw danych resekwencjonowania populacji podwojonych haploidów (DH) uzyskanych ze skrzyżowania dwóch linii kapusty chińskiej17. Następnie oceniliśmy lokalizację binmarkerów na mapach genetycznych, integrując je z odpowiednimi mapami fizycznymi. Z 892 binmarkerów w naszej asemblacji, 877 binmarkerów (98.3%) zostało zmapowanych na mapie genetycznej. Nasze złożenie zgadzało się z mapą genetyczną dla 801 binmarkerów (91.3%), co wskazuje na wysoką jakość wersji 3.0 (Rys.1; Tabela S7). Zauważyliśmy jednak, że 76 (8,7%) binmarkerów na chromosomach A05, A08 i A09 mapowało się do niejednoznacznych miejsc na mapie genetycznej. Regiony te zawierały powtarzające się sekwencje, zwłaszcza w regionach centromerycznych, jak opisano w poniższej analizie. Jednakże, te sprzeczne regiony zostały pokryte przez odczyty PacBio i/lub mapy BioNano; dane dla chromosomu A08 w v3.0 są pokazane jako przykład (Supplementary Figure S2).

Na mapie genetycznej v2.5 były 1092 binmarkery, a na mapie genetycznej v1.5 866 binmarkerów. Jednak udało nam się zmapować tylko 88,7% binmarkerów (969 z 1092) i 92,3% binmarkerów (799 z 866) na mapie genetycznej odpowiednio v2.5 i v1.5 (Supplementary Table S7). Stwierdziliśmy, że 15,1% binmarkerów (166 z 969) w v2.5 było rozbieżnych, w tym 146 binmarkerów z nieuporządkowanymi odległościami genetycznymi i fizycznymi w obrębie tego samego chromosomu (intrachromosom) i 20 binmarkerów z niespójnymi odległościami genetycznymi i fizycznymi na różnych chromosomach (interchromosom) (Supplementary Figure S3; Supplementary Table S7). Dla v1.5, 10,0% binmarkerów (80 z 799) było rozbieżnych, w tym 71 binmarkerów intrachromosomowych i 9 binmarkerów interchromosomowych (Supplementary Figure S4; Supplementary Table S7). Jednakże v3.0 zawierał najmniej sprzecznych binmarkerów intrachromosomalnych (8,7%, 76 z 877 binmarkerów) i żadnych rozbieżnych binmarkerów interchromosomalnych (Supplementary Table S7), wskazując, że genom B. rapa v3.0 ma wyższą przyległość niż dwa poprzednie złożenia. Łącznie, te niezależne walidacje sugerują, że B. rapa genome v3.0 ma najwyższą przyległość i najlepsze uporządkowanie rusztowań wśród trzech złożeń B. rapa.

Porównanie anotacji genomu

Przewidzieliśmy i anotowaliśmy modele genów, jak wcześniej opisano6. Zidentyfikowaliśmy łącznie 45 985 modeli genów kodujących białka w wersji 3.0, co stanowiło 14,74% zespołu genomu (Tabela 1). W naszym zespole 98,75% (45,411 z 45,985) genów zostało zaanotowanych na chromosomach, a tylko 1,25% (574 z 45,985) znajdowało się na rusztowaniach. Geny anotowane de novo w v3.0 zostały nazwane zgodnie ze standardem nomenklatury modeli genowych dla genomów referencyjnych Brassica (http://www.brassica.info/info/genome_annotation.php). Liczba modeli genów w nowym montażu jest wyższa niż w v1.5 (41,020 genów), ale niższa niż w v2.5 (48,826 genów) (Tabela 1). Aby dokładniej ocenić jakość anotacji, porównano ją z anotacją poprzednich złożeń przy użyciu BUSCO18, które oparte jest na benchmarku 1440 konserwatywnych genów roślinnych. Zidentyfikowano około 97,7% tych konserwatywnych genów roślinnych, a 1,7% wykryto jako fragmenty przedstawione w wersji 3.0 (Tabela UzupełniającaS11).

Analizę syntenii genomu przeprowadzono wśród trzech zespołów przy użyciu SynOrths19 w celu identyfikacji syntenicznych par genów i tandemowych tablic genów. W wersji 3.0 zidentyfikowano w sumie 2077 tandemowych tablic (odpowiadających 4963 tandemowym genom). Taką samą liczbę tandemowych tablic (2077 tablic odpowiadających 5004 genom) wykryto również w wersji v1.5. Ocena syntenii genomowej wykazała, że 1539 tandemowych tablic (odpowiadających 3757 genom) w v3.0 było syntenicznych z 1494 tandemowymi tablicami (odpowiadającymi 3670 genom) w v1.5. Jednakże w v2.5 zidentyfikowano więcej tablic tandemowych (3535 tablic, 8002 genów) (Tabela 1). Wykryliśmy luki w regionach zbędnych genów tandemowych w v2.5, podczas gdy w 3.0 i v1.5 nie znaleziono żadnych luk (Rys. 2a). Luki te mogą być wynikiem błędów w montażu powstałych w wyniku zamykania luk przy użyciu odczytów PacBio w v2.5, co z kolei doprowadziło do nieprawidłowej anotacji genów tandemowych. Dla innych genów tandemowych bez luk zaobserwowaliśmy, że pojedyncze geny w v3.0 i v1.5 były anotowane jako dwa lub więcej genów w v2.5 (Rys. 2b).

Po potraktowaniu każdej tablicy tandemowej jako pojedynczego locus genowego, pozostało 43 099 genów w v3.0, 44 359 genów w v2.5 i 38 093 genów w v1.5 (Tabela 1). Następnie przeprowadziliśmy analizę syntezy genów, która ujawniła, że 39 858 genów (92,48%) w v3.0 służyło jako odpowiedniki 40 442 (91,17%) i 35 464 genów (93,10%) odpowiednio w v2.5 i v1.5. Po porównaniu anotowanych genów z genami z wczesnych wersji, zidentyfikowaliśmy 3241 genów specyficznych dla wersji v3.0 w porównaniu do wersji v2.5 i v1.5. Spośród nich, 2380 genów było poparte dowodami z dopasowanych odczytów mRNA B. rapa (z BRAD, http://brassicadb.org/), a 2295 genów było poparte sekwencjami białek innych gatunków Brassicaceae (Tabela S12). W sumie 89,10% (2888 z 3214) genów specyficznych dla danej wersji w v3.0 było wspieranych przez dane mRNA B. rapa lub sekwencje białkowe innych gatunków Brassicaceae, podczas gdy tylko 10,90% (326 z 3214) genów nie było wspieranych.

Nowe zdarzenie ekspansji LTR-RT zidentyfikowane w zaktualizowanej asemblacji

Zanotowaliśmy TE w v3.0 przy użyciu tych samych metod, co poprzednio zgłoszone20. W sumie zidentyfikowano 235 683 TE z 1244 rodzin w v3.0, a 562 unikalne rodziny TE zostały znalezione w porównaniu z v2.5 i v1.5. W v3.0, TE stanowiły 37.51% (134 Mb) złożonego genomu, co było wyższą wartością niż w poprzednich złożeniach (32.30%, 126 Mb, v2.5; 25.44%, 72 Mb, v1.5)2,6. W naszej nowej asemblacji najliczniej występującymi TE są LTR-RT, których całkowita długość wynosi 57,64 Mb i które stanowią 16,32% zmontowanego genomu. Powtórzenia inne niż LTR-RT (LINE i SINE) stanowią 3,10% naszej asemblacji (Supplementary Figure S5). Wykryliśmy transpozony DNA odpowiadające 26,35 Mb, które stanowią 7,46% zmontowanego genomu (Supplementary Figure S5). Pełna lista zidentyfikowanych TE i powtórzeń w wersji 3.0 znajduje się w Tabeli S13. Ponadto zidentyfikowaliśmy 1231 miRNA, 1281 tRNA, 2865 rRNA i 3737 snRNA w genomie B. rapa v3.0 (Tabela uzupełniająca S19).

W naszym aktualnym złożeniu zanotowaliśmy więcej LTR-RT (57 Mb) w porównaniu do v2.5 (44 Mb) i v1.5 (18 Mb). Zidentyfikowaliśmy 51,062 nieuszkodzonych LTR-RT w v3.0. Dalsza analiza ujawniła, że 65,27% (33 672 z 51 602) nienaruszonych LTR-ów znajdowało się na dziesięciu chromosomach, podczas gdy 34,73% (17 922 z 51 602) nienaruszonych LTR-RT-ów znajdowało się na rusztowaniach niekotwiczonych. Używając tej samej metody6 , w wersji 3.0 zanotowano łącznie 13 318 nienaruszonych LTR-RT. Jednakże, w v2.5 i v1.5 było tylko 4129 i 801 nienaruszonych LTR-RT, odpowiednio6. Dalsza analiza ujawniła, że tylko 18,19% nienaruszonych LTR-RT (2423 z 13 318) znajdowało się na dziesięciu chromosomach, podczas gdy większość (81,81%, 10 895 z 13 318) nienaruszonych LTR-RT znajdowało się na rusztowaniach nie zakotwiczonych w v3.0.Czas insercji nienaruszonych LTR-RT obliczono zgodnie z wcześniejszym opisem4, co wskazuje, że genom B. rapa przeszedł trzy fale ekspansji LTR-RT od czasu odróżnienia się od B. oleracea (Rys. 3). Te nienaruszone LTR-RT miały średni wiek insercji 1,88 miliona lat temu (MYA), z medianą wieku insercji wynoszącą 1,59 MYA. Ponadto znaleźliśmy więcej nienaruszonych LTR-RT o różnych długościach w v3.0 w porównaniu z v2.5 i v1.5 (Supplementary Figure S6).

Dzięki tym nienaruszonym LTR-RT, zidentyfikowano nowe wydarzenie ekspansji LTR-RT w genomie B. rapa. Oznaczyliśmy 3155 nienaruszonych zdarzeń insercji LTR-RT od 0 MYA do 0,4 MYA jako „młodą ekspansję” o średniej długości 8135 bp i średniej dacie wstawienia 0,2 MYA; 2283 nienaruszonych zdarzeń insercji LTR-RT od 1,0 MYA do 1,4 MYA jako „średnią ekspansję” o średniej długości 11 902 bp i średniej dacie wstawienia 1,2 MYA; i 1444 nienaruszonych zdarzeń insercji LTR-RT od 3.0 MYA do 3,4 MYA jako „starożytna ekspansja” ze średnią długością 9823 bp i datą wstawienia (Rys. 3).Młoda i starożytna ekspansja ściśle odpowiadają wcześniej zidentyfikowanym zdarzeniom ekspansji; średnia ekspansja została po raz pierwszy zidentyfikowana w genomie B. rapa i ma podobny czas wstawienia jak nienaruszona ekspansja LTR-RT w B. oleracea. Ponadto, w wersji 3.0 zidentyfikowano 1778 LTR-RT podobnych do Ty1/Copia i 4179 LTR-RT podobnych do Ty3/Gypsy, czyli znacznie więcej niż w poprzednich asemblacjach (353 Ty1/Copia i 632 Ty3/Gypsy w wersji 2.5, 260 Ty1/Copia i 162 Ty3/Gypsy w wersji 1.5) (Tabela uzupełniająca S20; Rysunek uzupełniający S7, S8). Ogólnie rzecz biorąc, było więcej LTR-RT podobnych do Ty3/Gypsy niż LTR-RT podobnych do Ty1/Copia (Supplementary Table S20). W porównaniu z v2.5 i v1.5, Ty3/Gypsy-like LTR-RTs w v3.0 były wyraźnie zwiększone od 5 MYA (Supplementary Figure S7), podczas gdy Ty1/Copia-like LTR-RTs były zwiększone od 2,2 MYA (Supplementary Figure S8). Z drzew filogenetycznych stwierdziliśmy, że każda grupa LTR-RTs miała więcej kopii w v3.0 niż w v2.5 i v1.5 (Tabela uzupełniająca S21, S22; Rysunek uzupełniający S9, S10).

Bloki genomowe i centromery w genomie B. rapa

Badaliśmy relacje bloków genomowych używając zaktualizowanego złożenia v3.0. Aby zdefiniować bloki genomowe i centromery w genomie B. rapa v3.0, najpierw skonstruowaliśmy trzy subgenomy (LF, MF1 i MF2) w oparciu o relacje synteniczne między v3.0 i A. thaliana (Supplementary Figure S11; Supplementary Table S14). Wykryliśmy 71 z 72 (3 × 24) oczekiwanych bloków genomowych w v3.0, a większość z nich była ułożona zgodnie z blokami podanymi wcześniej w ref.10 (Ryc. 4; Tabela uzupełniająca S15). W v3.0, dwa nowe pofragmentowane bloki genomowe F (LF) i F (MF1) zostały zidentyfikowane na chromosomach A01 i A05 i nie były obserwowane w ref.10. Nie udało nam się wykryć dwóch wcześniej opisanych, bardzo małych bloków genomowych w v3.0, bloku C (MF2) na chromosomie A07 i bloku B (MF1) na chromosomie A08 (ref.10). Jednak w naszym montażu bloki genomowe N/M (MF1), O/P (LF) i A/C(LF) były ułożone odpowiednio na chromosomach A01, A09 i A10, podczas gdy w ref.10 były one ułożone po przeciwnych stronach. Trzy małe sąsiadujące bloki genomowe (S (MF2), T (MF2) i B (MF1)) na chromosomie A08 w wersji 3.0 były uporządkowane jako S/T/B, podczas gdy w ref.10 były one ułożone jako T/B/S..

Porównaliśmy również orientację bloków genomowych w v3.0 z tą w ref.10. Bloki genomu W (MF2) na chromosomie A02, jak również G (LF) i E (LF) na chromosomie A07, okazały się być odwrócone w stosunku do innych bloków, które pochodzą z pojedynczego chromosomu ACK. Jednakże orientacja bloku genomu P (LF) na chromosomie A09 i trzech bloków V w v3.0 była w kierunku do przodu, podczas gdy były one odwrócone w ref.10. Wyniki te zostały dodatkowo poparte przez mapy genetyczne odpowiednio v3.0 i v1.5.

Dokładnie określiliśmy położenie centromerów wszystkich chromosomów w v3.0. Przesiewając wcześniej określone sekwencje powtórzeń centromerowych, w tym centromerowe powtórzenia satelitarne CentBr, CRB, TR238 i PCRBr21,22,23, zidentyfikowaliśmy sygnały dla wszystkich 21 paleocentromerowych regionów w v3.0, podczas gdy trzy paleocentromerowe regiony nie zostały wykryte w ref.10 (ryc. 4, tabela uzupełniająca S16). Analiza paleocentromerowa wykazała, że wszystkie dziesięć istniejących centromerów B. rapa zostało odziedziczonych z 21 paleocentromerów. W wersji 3.0, centromery chromosomów A01, A03, A04, A05, A06, A07 i A10 miały te same bloki genomu otaczające odpowiadające im centromery, jak podano w ref.10 (Ryc. 4). Jednakże centromer na chromosomie A02 znajdował się pomiędzy blokami genomowymi P (MF2) i V (MF1), a centromer na chromosomie A09 znajdował się pomiędzy blokami genomowymi P (LF) i B (LF), podczas gdy w ref.10 były one uznane za paleocentromery (Rys. 4). Centromer na chromosomie A08 znajdował się pomiędzy blokami genomowymi T (MF2) i B (MF1), a nie pomiędzy blokami genomowymi C (MF1) i T (MF2), jak podano w ref.10. Ponadto w v3.0 wykryto 1188 genów w regionach centromerowych, podczas gdy w ref.10 wykryto tylko 740 genów (Tabela uzupełniająca S17).

Aby ocenić nasze złożenie w odniesieniu do centromerów w v3.0, przeanalizowaliśmy cechy sekwencji regionów centromerowych. Stwierdziliśmy, że znacznie większa liczba TE i powtórzeń specyficznych dla centromeru została zmapowana do regionów centromerowych niż do innych części chromosomów, a gęstość genów i wskaźnik rekombinacji były wyraźnie niższe w regionach centromerowych opisanych w wersji 3.0 (ryc. 5). Ponadto, w regionach centromerowych w v3.0 wykryto więcej powtórzeń specyficznych dla centromeru w porównaniu z tymi opisanymi w ref.10 (Supplementary Table S17).

.