Assembly and annotation of C. kanehirae
SCT jest diploidalny (2n = 24; Supplementary Fig. 1a) z szacowanym rozmiarem genomu 823,7 ± 58,2 Mb/1 C (Supplementary Fig. 1b i 2). Stworzyliśmy asemblację pochodzącą wyłącznie z 85× PacBio long reads (read N50 = 11,1 kb; contig N50 = 0,9 Mb) obejmującą 728,3 Mb. Sekwencje konsensusowe tej asemblacji zostały skorygowane przy użyciu 141× odczytów Illumina i dalej rusztowane przy użyciu 207× odtworzonej chromatyny 'Chicago’ i 204× sparowanych odczytów Hi-C przy użyciu potoku HiRise (Supplementary Fig. 3). Uzyskano ostateczne, zintegrowane złożenie 730,7 Mb w 2153 rusztowaniach, obejmujące 91,3% szacowanego rozmiaru genomu z cytometrii przepływowej. Ostateczne rusztowanie N50 wynosiło 50,4 Mb, z ponad 90% w 12 pseudomolekułach (Tabela uzupełniająca 1), przypuszczalnie odpowiadających 12 chromosomom SCT.
Używając kombinacji wsparcia homologii białek roślin referencyjnych i sekwencjonowania transkryptomu pochodzącego z różnych tkanek (Supplementary Fig. 1c i Tabela 2) oraz przewidywania genów ab initio, 27 899 modeli genów kodujących białka zostało przypisanych przy użyciu potoku MAKER218 (Tabela uzupełniająca 1). Spośród nich 93,7% okazało się homologiczne z białkami w bazie TrEMBL, a 50% można było przypisać do terminów Gene Ontology przy użyciu eggNOG-mapper19. Proteom został oszacowany jako co najmniej 89% kompletny na podstawie oceny BUSCO20 (benchmarking universal single-copy orthologs), co jest porównywalne z innymi sekwencjonowanymi gatunkami roślin (Supplementary Table 1). Klasteryzacja Orthofinder21 modeli genów SCT z modelami pochodzącymi z 12 różnych genomów roślin nasiennych dała 20 658 grup ortologicznych (Tabela Uzupełniająca 3). 24 148 genów SCT (86,56%) należało do grup ortologicznych z ortologami z co najmniej jednego innego gatunku roślin. 3 744 modele genów nie były ortologiczne z innymi, a tylko 210 genów należało do 48 grup ortologicznych specyficznych dla SCT. W sumie sugerują one, że zróżnicowanie fenotypowe u magnoliowatych może być napędzane przez narodziny de novo genów specyficznych dla gatunku i ekspansję istniejących rodzin genów.
Charakterystyka genomu
Zidentyfikowaliśmy 3 950 027 heterozygotycznych miejsc biallelicznych w genomie SCT, co odpowiada średniej heterozygotyczności 0,54% (jeden heterozygotyczny polimorfizm pojedynczego nukleotydu (SNP) na 185 bp). Częstości alternatywnych (niereferencyjnych) alleli w tych miejscach miały główny szczyt około 50%, co jest zgodne z faktem, że SCT jest diploidalny i nie ma dowodów na niedawną aneuploidię (Supplementary Fig. 4). Przestrzenne rozmieszczenie miejsc heterozygotycznych było bardzo zróżnicowane, przy czym 23,9% genomu wykazywało mniej niż 1 locus SNP na kb w porównaniu do 10% genomu z co najmniej 12,6 loci SNP na kb. Przebiegi regionów homozygotyczności wydawały się być rozmieszczone losowo na chromosomach SCT, osiągając maksimum 20,2 Mb w rusztowaniu 11 (Ryc. 1a). Takie długie przebiegi regionów homozygotyczności mają takie samo pokrycie sekwencji jak reszta genomu (Supplementary Fig. 5) i mogą być związane z selektywnym wymazywaniem, chowem wsobnym lub niedawnym wąskim gardłem populacji. Geny zlokalizowane w tych regionach homozygotyczności okazały się być wzbogacone w proces biosyntezy ligniny i metabolizm galaktozy (Tabela 4), co sugeruje ich potencjalną rolę w tworzeniu kompleksów ligninowo-węglowodanowych22. Analiza parami sekwencyjnie markowskiej koalescencji23 (PSMC) oparta na gęstości heterozygotycznych SNP wskazała na ciągłe zmniejszanie się efektywnej wielkości populacji w ciągu ostatnich 9 milionów lat (ryc. 1b), z możliwym wąskim gardłem zbiegającym się z przesunięciem klimatycznym w połowie plejstocenu 0,9 miliona lat temu (Ma). Takie wzorce mogą odzwierciedlać złożoną historię populacji SCT związaną z geologiczną historią Tajwanu, w tym wypiętrzanie i formowanie wyspy w późnym miocenie (9 Ma), po którym nastąpiła budowa gór 5-6 Ma, odpowiednio24.
a, Liczba heterozygotycznych biallelicznych SNPs na 100-kb nienakładających się okien jest wykreślona wzdłuż największych 12 rusztowań. Wykluczono indeksy. b, Historia efektywnej wielkości populacji została wyznaczona przy użyciu metody PSMC. Wykonano sto bootstrapów, a marginesy zaznaczono kolorem jasnoczerwonym. c, Dla każdego nienakładającego się okna 100-kb przedstawiono rozkład od góry do dołu: gęstość genów (procent nukleotydów z przewidywanym modelem), transkryptom (procent nukleotydów z dowodami mapowania transkryptomu) i trzy różne klasy sekwencji powtarzalnych (procent nukleotydów z adnotacją elementu transpozycyjnego). Czerwona litera T oznacza obecność telomerycznego klastra powtórzeń na końcu rusztowania. LINE, long interspersed nuclear element.
Części transpozycyjne i powtórzenia interspersyjne stanowiły 48% zespołu genomu (Tabela 5). Większość elementów transpozycyjnych należała do retrotranspozonów o długim terminalnym powtórzeniu (LTR) (25,53%), a następnie do elementów transpozycyjnych DNA (12,67%). Wśród LTR-ów, 40,75% i 23,88% retrotranspozonów należało odpowiednio do Ty3/Gypsy i Ty1/Copia (Tabela 5). Filogeneza domeny odwrotnej transkryptazy wykazała, że większość kopii Ty3/Gypsy tworzyła odrębny klad (20 092 kopii), przypuszczalnie w wyniku niedawnej ekspansji i proliferacji, podczas gdy elementy Ty1/Copia zgrupowane były w dwóch siostrzanych kladach (7 229 i 2 950 kopii) (Supplementary Fig. 6). Z wyjątkiem dwóch rusztowań, zarówno Ty3/Gypsy, jak i Ty1/Copia elementy transpozycyjne LTR były skupione w centrach pericentromerycznych 12 największych rusztowań (ryc. 1c i ryc. 7). Ponadto, regiony wzbogacone w LTR (zdefiniowane przez 100 kb z ponad 50% udziałem elementów transpozycyjnych klasy LTR) miały średnio o 35% większe pokrycie niż reszta genomu (Ryc. 1c i Suplementary Fig. 8), co sugeruje, że powtórzenia te zostały zawalone w montażu i mogły przyczynić się do różnic w szacunkach wielkości genomu na podstawie cytometrii przepływowej i k-mer. Zawartość sekwencji kodujących w SCT jest podobna do innych genomów okrytozalążkowych uwzględnionych w naszych analizach (Tabela 1), podczas gdy introny są nieco dłuższe w SCT z powodu większej gęstości elementów transpozycyjnych (P < 0,001, test Wilcoxon rank-sum; Supplementary Fig. 9).
Jak opisano dla innych genomów roślinnych25, rusztowania na poziomie chromosomów SCT wykazują niską gęstość genów kodujących białka i wysoką gęstość elementów transpozycyjnych w centrach chromosomów, a także zwiększoną gęstość genów w kierunku końców chromosomów (Fig. 1c). Zidentyfikowaliśmy skupiska przypuszczalnego heptameru subtelomerowego, TTTAGGG, rozciągające się na długość aż 2547 kopii, co wskazuje na występowanie powtórzeń telomerowych u roślin26 (Tabela 6). Ponadto, odkryto 687 kb sekwencji podobnych do jądrowego plastydowego DNA (NUPTs), o średniej długości około 202,8 bp (Tabela 7). Wśród NUPT SCT przeważały krótkie fragmenty, przy czym 96% zidentyfikowanych NUPT miało długość mniejszą niż 500 bp (Tabela 8). Najdłuższy NUPT ma długość ~20 kb i jest synteniczny z 99,7% identycznością z fragmentem plastomu SCT, który zawiera siedem genów kodujących białka i pięć genów tRNA (Supplementary Fig. 10).
Phylogenomiczne umiejscowienie C. kanehirae siostrzanym do eudicots
Aby rozwiązać długotrwałą debatę na temat filogenetycznego umiejscowienia magnoliowatych w stosunku do innych głównych linii roślin kwitnących, skonstruowaliśmy drzewo filogenetyczne oparte na 211 ściśle jednokopijnych zestawach ortologów (to jest jeden i tylko jeden homolog we wszystkich gatunkach) zidentyfikowanych przez OrthoFinder21 dookoła rodziny genów wszystkich modeli genów z SCT i 12 innych genomów roślin nasiennych (patrz Metody). Pojedyncze drzewo gatunkowe zostało odtworzone poprzez analizę maksymalnego prawdopodobieństwa27 konkatenowanej supermacierzy dopasowań genów pojedynczej kopii oraz analizę opartą na koalescencji z wykorzystaniem drzew 211 genów28 (Ryc. 2 i Suplementary Fig. 11). SCT, reprezentujący linię magnoliowatych, został umieszczony jako siostrzany do kladu eudicot (Rys. 2). Topologia ta pozostała stabilna, gdy uwzględniliśmy zestaw danych transkryptomowych dodatkowych 22 gatunków magnoliowatych z inicjatywy 1000 roślin29 (1KP), choć uzyskano niższe wsparcie bootstrapowe (Supplementary Fig. 12). Używając MCMCtree30 z kalibracjami kopalnymi, obliczyliśmy 95% przedział ufności dla czasu dywergencji między magnoliowcami a eudikotami na 136,0-209,4 Ma (ryc. 2), który pokrywa się z dwoma innymi niedawnymi szacunkami (114,8-164,1 Ma31 i 118,9-149,9 Ma32).
Rozszerzenie i kurczenie się rodziny genów oznaczono odpowiednio liczbami obok znaków plus i minus. Zielone liczby w nawiasach oznaczają szacunkowe czasy dywergencji. Wsparcie bootstrapowe wszystkich węzłów wynosiło 100, chyba że podano inaczej.
Analiza syntenii/duplikacja całego genomu
Poprzednie badania danych znaczników sekwencji wyrażonych wnioskowały o genomowej duplikacji w obrębie magnoliowatych przed dywergencją Magnoliales i Laurales10, ale testowanie tej hipotezy w oparciu o syntenię nie było możliwe bez zmontowanego genomu magnoliowatych. W sumie zidentyfikowano 16 498 par genów w 992 blokach syntenicznych, stanowiących 72,7% złożenia genomu SCT. Spośród tych wewnątrzgenomowych bloków syntenicznych 72,3% było syntenicznych z więcej niż jednym miejscem w genomie, co sugeruje, że w rodowodzie SCT doszło do więcej niż jednej duplikacji całego genomu (WGD) (Ryc. 3a). Dwie rundy starożytnej WGD były sugerowane przez rozległą syntenię między parami regionów chromosomalnych i znaczące, ale mniej synteniczne sparowanie każdego regionu z dwoma dodatkowymi segmentami genomowymi (Supplementary Fig. 13). Bloki syntenii 12 największych rusztowań SCT przyporządkowano do pięciu klastrów, które mogą odpowiadać chromosomom przodków pre-WGD (Fig. 3a, Supplementary Fig. 13 i Supplementary Note).
a, Schematyczne przedstawienie wewnątrzgenomowych relacji między 637 blokami syntenii w genomie SCT. Bloki syntenii (oznaczone przez bloki brzoskwiniowe) zostały jednoznacznie przypisane do pięciu klastrów powiązań reprezentujących starożytne kariotypy i są oznaczone kolorem. Fioletowe bloki oznaczają blok syntenii przypisany w pierwszej grupie sprzężeń (patrz także Supplementary Fig. 13). b, Schematyczne przedstawienie pierwszej grupy sprzężeń w genomie SCT i odpowiadających im relacji u A. trichopoda.
Amborella trichopoda jest jedynym gatunkiem reprezentującym linię siostrzaną do wszystkich innych wymarłych roślin okrytozalążkowych i nie ma dowodów na WGD od czasu dywergencji z ostatniego wspólnego przodka wymarłych linii roślin kwitnących33. Aby potwierdzić, że dwie rundy WGD miały miejsce w rodowodzie SCT po rozejściu się linii prowadzących do SCT i A. trichopoda, oceniliśmy syntenię między tymi dwoma genomami. Zgodnie z naszą hipotezą, od jednego do czterech segmentów genomu SCT było wyrównanych do jednego regionu w genomie A. trichopoda (Fig. 3b i Supplementary Fig. 14).
Aby dokładniej określić czas dwóch rund WGD widocznych w genomie SCT, oszacowano wewnątrzgatunkowe i międzygatunkowe rozkłady homologii Ks (substytucje synonimiczne na miejsce synonimiczne). Wewnątrzgenomowe duplikaty SCT wykazały dwa szczyty na poziomie 0,46 i 0,76 (Ryc. 4a), odpowiadające dwóm zdarzeniom WGD. Na podstawie tych dwóch szczytów mogliśmy wnioskować o ewolucji kariotypu, organizując zgrupowane bloki syntenii w cztery grupy, które prawdopodobnie pochodzą z jednego z pięciu chromosomów sprzed WGD (Supplementary Fig. 15). Porównanie między Aquilegia coerulea (Ranunculales, linia siostrzana do wszystkich innych wymarłych eudicots33) i ortologami SCT ujawniło wyraźny pik wokół Ks = 1,41 (Ryc. 4a), podczas gdy wewnątrzgenomowy duplikat Aquilegia wynosił około Ks = 1, sugerując niezależne WGD po dywergencji linii prowadzących do SCT i Aquilegia. Dostępność transkryptomu 17 Laurales + Magnoliales z 1KP29 pozwoliła nam sprawdzić hipotezę o czasie wystąpienia WGDs widocznych w genomie SCT8. Rozkład Ks dla pięciu z sześciu dostępnych gatunków z rodziny Lauraceae ujawnił dwa piki (Rys. 4b i Supplementary Fig. 16), podobnie jak w przypadku rozkładu Ks SCT (Rys. 4a) i odpowiadające dwóm opartym na syntenii wnioskom o WGD w rodowodzie SCT (Rys. 3 i Supplementary Fig. 15). U pozostałych gatunków Laurales i Magnoliales zaobserwowano tylko jeden pik Ks, co sugeruje, że w rodowodzie tych gatunków wystąpiło tylko jedno zdarzenie WGD (Supplementary Figs. 17 i 18). Szczyt Ks widoczny w danych Aquilegia jest prawdopodobnie przypisywany WGD w obrębie Ranunculales długo po dywergencji eudicots i magnoliids (Fig. 4a).
a, Duplikaty par ortologów zidentyfikowane w blokach syntenicznych w obrębie SCT, A. coerulea oraz między SCT i A. coerulea. b, Ks wewnątrzgenomowych duplikatów par w obrębie Lauraceae i Magnoliales w projekcie 1KP29. Linie przerywane oznaczają dwa szczyty Ks zaobserwowane w SCT. Brązowe i szare linie oznaczają odpowiednio rozkład Ks u SCT i innych Lauraceae.
Specjalizacja proteomu magnoliowatych
Dążymy do zidentyfikowania genów i domen białkowych specyficznych dla SCT poprzez anotację domen rodziny białek (Pfam) i ocenę ich rozmieszczenia w 13 genomach roślin nasiennych włączonych do naszych analiz filogenomicznych. Zgodnie z obserwacją, że istnieje bardzo niewiele grup ortologicznych specyficznych dla SCT, analiza składowych głównych zawartości domen Pfam połączyła SCT z monokotami i eudikotami, przy czym dwie pierwsze składowe główne oddzieliły od tej grupy gymnospermy i A. trichopoda (Supplementary Fig. 19a). Stwierdzono znaczne nakładanie się SCT, gatunków eudikota i monokota, co sugeruje znaczne zróżnicowanie funkcjonalne od czasu rozdzielenia się tych trzech linii. SCT wykazał również znaczące wzbogacenie i redukcję odpowiednio 111 i 34 domen białkowych w porównaniu z innymi gatunkami roślin (Supplementary Fig. 19b i Supplementary Table 9). Wzbogacenie domen białkowych obejmowało karboksy-końcową domenę syntazy terpenowej (TPS) zaangażowaną w reakcje obronne oraz powtórzenia bogate w leucynę (628 versus 334,4) w efektywność transpiracji roślin34. Co ciekawe, stwierdziliśmy, że SCT posiada 21 kopii czynnika transkrypcyjnego EIN3/EIN3-like (EIL), więcej niż poprzednio odnotowana maksymalna liczba 17 kopii w genomie banana (Musa acuminata)35. EIL inicjują sygnalizację etylenową poprzez aktywację czynnika odpowiedzi etylenowej (ERF), którego ekspansję stwierdziliśmy również w SCT (150 kopii w porównaniu do średniej 68,3 kopii z dziewięciu gatunków opisanych w ref. 35; Supplementary Fig. 20). ERF odpowiada i pozytywnie moduluje biosyntezę sygnałów fitohormonalnych, w tym etylenu36. Ekspresja ERF została włączona w pozytywną modulację rozwoju roślin od dojrzewania owoców35 do wzrostu wtórnego w tworzeniu drewna37, jak również w zwiększoną odporność na czynniki abiotyczne38 lub biotyczne39. Tak więc, ekspansja EILs w SCT może stymulować ERF, prowadząc do różnych regulacji efektorów downstream, które skutkują cechami specyficznymi dla SCT.
Następnie oceniliśmy ekspansje i kurczenie się grup ortologicznych w filogenezie roślin nasiennych (Ryc. 2). Ewolucja wielkości rodziny genów była dynamiczna w całej filogenezie, a gałąź prowadząca do SCT nie wykazywała znacząco różnej liczby ekspansji i kurczenia się. Wzbogacenie terminów z zakresu ontologii genów ujawniło albo różne rodziny genów pełniące wspólne funkcje, albo rodziny pojedynczych genów ulegające znacznej ekspansji (tabele uzupełniające 10 i 11). Na przykład, rozszerzeni członkowie genów odporności roślin (R) sumują się do „odpowiedzi nadwrażliwości typu roślinnego” (Tabela uzupełniająca 10). Natomiast wzbogacone terminy Gene Ontology z zakontraktowanych rodzin genów gałęzi SCT (Tabela uzupełniająca 11) zawierają członków transporterów ABC, syntetazy kwasu indolilo-3-octowego-amido, endotransglukozylazy/hydrolazy ksyloglukanu i białka reagującego na auksynę, z których wszystkie są częścią „odpowiedzi na auksynę”.
Geny R
Anotacja genomu SCT zawierała 387 modeli genów R, z których 82% należy do typów nucleotide-binding site leucine-rich repeat (NBS-LRR) lub coiled-coil NBS-LRR. Wynik ten jest zgodny z wcześniejszym doniesieniem, że LRR jest jedną z najliczniej występujących domen białkowych u roślin i jest wysoce prawdopodobne, że SCT jest w stanie rozpoznawać i zwalczać produkty genów awirulencji (Avr) patogenów40. Wśród badanych 13 genomów, SCT posiada największą liczbę genów R wśród roślin nieuprawnych (Supplementary Fig. 21). Drzewo filogenetyczne zbudowane z 2465 domen NBS sugeruje również, że klady w obrębie tej rodziny genów różnicowały się niezależnie w obrębie eudikotów, monokotów i magnoliowatych. Co ciekawe, najbardziej zróżnicowane klady genów NBS SCT były siostrzane wobec zdepauperyzowanych kladów genów NBS eudicota (Supplementary Fig. 22).
Rodzina genów TPS
Jedną z najbardziej uderzających cech genomu SCT jest duża liczba genów TPS (CkTPS). Łącznie przewidziano i zanotowano 101 genów CkTPS, co stanowi największą liczbę w jakimkolwiek innym genomie. Dzięki włączeniu zestawu danych transkryptomu dwóch kolejnych gatunków magnoliowatych (P. americana i Saruma henryi), analizy filogenetyczne TPS z 15 gatunków umieściły geny CkTPS wśród sześciu z siedmiu podrodzin genów TPS, które zostały opisane dla roślin nasiennych41 (Ryc. 5, Tabela 1 i Ryc. 23-28). Geny CkTPS umieszczone w podrodzinach TPS-c (2) i TPS-e (5) prawdopodobnie kodują syntazy diterpenowe, takie jak syntaza difosforanu kopalilu i syntaza ent-kaurenu42. Są to kluczowe enzymy katalizujące powstawanie 20-węglowych izoprenoidów (określanych zbiorczo jako diterpenoidy; C20s), o których sądzono, że są specyficzne dla eudikoty41 i pełnią podstawowe funkcje, takie jak regulacja pierwotnego metabolizmu roślin. Pozostałe 94 przewidywane geny CkTPS prawdopodobnie kodują 10-węglowe syntazy monoterpenowe (C10), 15-węglowe syntazy seskwiterpenowe (C15) i dodatkowe 20-węglowe syntazy diterpenowe (C20) (Tabela 1). Z odpowiednio 25 i 58 homologami, podrodziny TPS-a i TPS-b są najbardziej zróżnicowane w SCT, przypuszczalnie przyczyniając się do masowej i mieszanej produkcji lotnych C15s i C10s43.
Drzewo filogenetyczne zostało skonstruowane przy użyciu putatywnych lub scharakteryzowanych genów TPS z 13 sekwencjonowanych genomów roślin lądowych i dwóch magnoliowatych z dostępnymi danymi transkryptomicznymi.
Na uwagę zasługuje fakt, że drzewo genów TPS rozwiązało specyficzne dla Laurazji klady genów TPS w obrębie podrodzin genów TPS-a, TPS-b, TPS-f i TPS-g (ryc. uzupełniające 23, 24, 27 i 28). Ten wzór duplikacji genów TPS u wspólnego przodka Persea i Cinnamomum, a następnie ich retencji może wskazywać na subfunkcjonalizację lub neofunkcjonalizację zduplikowanych genów TPS w obrębie Lauraceae. W analizach zidentyfikowano także specyficzną dla magnoliowatych podkladę w podrodzinie TPS-a, obejmującą więcej magnoliowatych genów TPS o scharakteryzowanych funkcjach (Supplementary Fig. 23). Rzeczywiście, wykryliśmy pozytywną selekcję w specyficznych dla Laurazji podkladach TPS-f -I i -II, co sugeruje funkcjonalną dywergencję (Supplementary Table 13). Łącznie dane te wskazują na rosnące zróżnicowanie genów TPS magnoliowatych zarówno przed, jak i po powstaniu rodziny Lauraceae.
CkGeny TPS nie są równomiernie rozmieszczone na chromosomach (Tabela uzupełniająca 12), a skupiska członków poszczególnych podrodzin zaobserwowano jako tandemowe duplikaty (Ryc. uzupełniająca 29). Siedemdziesiąt sześć genów TPS zaobserwowano w największych 12 rusztowaniach SCT. Spośród nich, 60,5% (46 kopii) należących do różnych podrodzin znaleziono w regionie 0,5-15 Mb i 22,0-24,5 Mb odpowiednio w rusztowaniach 7 i 10 (Supplementary Fig. 29). Rusztowanie 7 zawiera 29 genów CkTPS należących do kilku podrodzin, w tym wszystkie z ośmiu CkTPS-a, 12 CkTPS-b, pięciu CkTPS-e i trzech CkTPS-f (ryc. 29). Natomiast tylko dwóch członków CkTPS-c rezyduje w rusztowaniu 1. Dwadzieścia cztery geny CkTPS zlokalizowane są w innych, mniejszych rusztowaniach, z których 22 kodują podrodzinę TPS-b (ryc. 24). Niektóre z tych podrodzin, zlokalizowane na rusztowaniach 7 i 10, znajdują się w fizycznym sąsiedztwie (ryc. 29). Na przykład 3 z 11 członków podrodziny TPS-b-Lau III sąsiadowały z 4 z 11 członków podrodziny TPS-b-Lau V (Supplementary Fig. 29), podczas gdy inni członkowie podrodziny znajdowali się nie w odpowiadających im regionach syntenicznych, lecz w innych miejscach genomu (Supplementary Fig. 30). Geny należące do tego klastra nie były zgrupowane razem w odpowiadającej im filogenezie podrodziny (Supplementary Fig. 30), co sugeruje, że ich ułożenie mogło nastąpić później niż ostatnie zdarzenie WGD.