Genome sequencing, assembly and annotation

Sekwencjonowaliśmy O. sinensis z dystryktu Nyingchi w Tybecie, Chiny. Przeprowadziliśmy analizę WGS za pomocą platform sekwencjonowania następnej generacji Roche 454 i Illumina HiSeq 2000. Wygenerowaliśmy zestawy czystych danych sekwencyjnych o wielkości ~5,4 Gb, co dało odpowiednio około 45,1-krotne pokrycie genomu (Tabela S1). Oszacowaliśmy, że rozmiar genomu wynosi ~ 124,08 Mb i ~ 119,8 Mb w oparciu o cytometrię przepływową i 17-merowy rozkład głębokości zsekwencjonowanych odczytów, odpowiednio (Supplemental Figures S1-2 i Supplemental Table S2). Genom O. sinensis został najpierw poskładany z długich odczytów Roche 454 przy użyciu Newblera16, a następnie posortowano wstępnie poskładane kontigi z odczytami sekwencjonowania parami par Illumina przy użyciu SSPACE17. Ostatecznie dało to złożenie genomu o wielkości ~116,4 Mb, które obejmuje ~97% szacowanego rozmiaru genomu i zawiera 156 rusztowań (>2 Kb) o wartości ScafN50 wynoszącej ~3 Mb i 9,141 kontigów (N50 = 21,423 bp) (Tabela 1 i Tabele uzupełniające S3-4). Aby zweryfikować jakość złożenia genomu, najpierw wyrównaliśmy wszystkie DNA i znaczniki sekwencji wyrażonych (ESTs) O. sinensis dostępne w publicznych bazach danych i uzyskaliśmy wskaźniki mapowania odpowiednio 98,85% i 95,33% (Tabela uzupełniająca S5). Po drugie, zmapowaliśmy wszystkie czyste odczyty Roche 454 (~1,84 Gb) do zmontowanych sekwencji genomu, wykazując niemal idealne dopasowanie ze wskaźnikiem mapowania 99,01% (Tabela S5). Po trzecie, zmontowane przez nas transkrypty wykazały dobre dopasowanie do zmontowanego genomu; z 11 742 transkryptów 91,29% zostało zmapowanych (pokrycie transkryptów ≥80% i identyczność ≥90%; Tabela S5). Wreszcie, oceniliśmy kompletność naszego złożenia O. sinensis przy użyciu BUSCO18; 94,0%, 4,0% i 1,8% z 1315 oczekiwanych konserwatywnych genów Ascomycota BUSCO zostało zidentyfikowanych odpowiednio jako kompletne, fragmentaryczne i brakujące w naszym złożeniu O. sinensis (Supplemental Table S5).

W celu wspomagania przewidywania genów wygenerowaliśmy ~15.05 Gb danych sekwencjonowania RNA (RNA-Seq) uzyskanych z łącznie sześciu bibliotek reprezentujących trzy główne stadia rozwojowe (Supplemental Figure S4 i Supplemental Tables S6,7). W połączeniu z predykcją ab initio, wyrównaniami białek i EST, łączeniem EvidenceModeler i dalszym filtrowaniem, zdefiniowaliśmy 7,939 genów kodujących białka (Tabela 1 i Uzupełniająca Tabela S8). Spośród tych przewidywanych genów, około 97,0% i 71,51% mogło być funkcjonalnie sklasyfikowanych i popartych danymi RNA-Seq, odpowiednio (Tabele uzupełniające S9-11). Używając BUSCO z linii Ascomycota, stwierdziliśmy, że 94,4%, 3,6%, 1,8% i 0,2% genów było odpowiednio kompletnych, fragmentarycznych, brakujących i zduplikowanych, co wskazuje na dobrą jakość naszej anotacji genów (Tabela S11). Przeprowadziliśmy również wyszukiwanie homologów i anotowaliśmy geny niekodującego RNA (ncRNA), uzyskując 146 genów transferowego RNA (tRNA), 33 geny rybosomalnego RNA (rRNA), 70 genów małego nukleolarnego RNA (snoRNA) i 15 genów małego jądrowego RNA (snRNA) (Supplemental Figure S6 i Supplemental Table S12). Anotacja sekwencji powtórzonych wykazała, że elementy transpozycyjne (TE) stanowiły około 74,67% zmontowanego genomu i 80,07% surowych odczytów, co wskazuje, że ~5,45% niezmontowanego genomu składa się z TE (Supplemental Tables S13-14). Zawartość GC wynosiła 43,09% w całym genomie i 61,49% w sekwencjach kodujących (Supplemental Figure S3; Supplemental Tables S4 i S8). Zanotowaliśmy 8,918 powtórzeń sekwencji prostej, które będą stanowiły cenne markery genetyczne pomocne w przyszłych programach hodowlanych chińskiego gąsienicowego grzyba (Supplemental Tables S15-16 i Supplemental Figure S7).

Retrotransposon-driven genome expansion and massive removal of non-collinear genes

Porównanie rozmiarów genomu wykazało, że genom O. sinensis był prawie 3,4-krotnie większy niż innych entomopatogenicznych grzybów z rodziny Hypocreales (Supplemental Table S17 i Supplemental Figure S8A). Analiza sekwencji powtórzeń wykazała, że ekspansja ta była spowodowana przede wszystkim szybkim namnażaniem się elementów transpozycyjnych. Około 74,67% zespołu genomu O. sinensis składało się z sekwencji powtarzalnych (Tabele uzupełniające S13-14), wyjątkowo większych niż te zgłaszane u Metarhizium anisopliae (~0,98%)19, Metarhizium acridum (~1,52%)19, Cordyceps militaris (~3,04%)20 i Beauveria bassiana (~2,03%)21 (P < 4,822e-07) (Supplemental Figure S8B). Elementy MULE były szczególnie obfite, stanowiąc ~ 1,6% (~ 1,9 Mb) genomu O. sinensis i ponad 59% transpozonów DNA w O. sinensis. Retrotranspozony, głównie retrotranspozony o długim okresie powtarzania (LTR), stanowiły ~59,76% genomu O. sinensis, a ich proliferacja na dużą skalę nastąpiła około ~38 milionów lat temu (MYA) (Supplemental Figure S9).

W przeciwieństwie do szybkiej amplifikacji retrotranspozonów LTR napędzających ekspansję genomu O. sinensis, inną niezwykłą cechą jest dramatyczna utrata genów kodujących białka w linii O. sinensis w porównaniu z innymi grzybami entomopatogenicznymi. W porównaniu z całkowitą liczbą 7 939 genów kodujących białka u O. sinensis, u innych grzybów entomopatogenicznych, np. Metarhizium anisopliae (10 582)19, Metarhizium acridum (9 849)19, Cordyceps militaris (9 684)20, Beauveria bassiana (10 366)21 i Tolypocladium inflatum (9 998)22, było średnio ponad 10 095 genów (tab. 1). O takiej redukcji liczby genów świadczyła również identyfikacja genów niekolinearnych oraz analiza porównawcza bloków syntenicznych pomiędzy genomami O. sinensis i C. militaris. Zidentyfikowaliśmy łącznie 308 bloków syntenicznych, które obejmują prawie 72,7% (~23,4 Mb w C. militaris vs. ~43,5 Mb w O. sinensis) genomu C. militaris (Ryc. 1A; Suplemental Figure S10 i Supplemental Tables S18-19). Spośród tych syntenicznych regionów genomowych odnotowano spadek liczby genów niekolinearnych u O. sinensis (2 127) w porównaniu z C. militaris (3 259), ale wzrost liczby sekwencji powtórzonych (23,8 Mb u O. sinensis vs. 0,40 Mb u C. militaris) (Ryc. 1B i Uzupełniająca Tabela S19). Funkcjonalna anotacja 2,468 genów, które zostały utracone w O. sinensis wykazała, że były one głównie zaangażowane w metabolizm aminokwasów, takich jak biosynteza aminokwasów (ko01230), metabolizm argininy i proliny (ko00330) oraz metabolizm tyrozyny (ko00350) (Supplemental Figure S11 i Supplemental Table S20). Warto zauważyć, że prawie 81% sekwencji powtórzonych w tych 308 syntenicznych blokach było retrotranspozonami LTR, z których 40,4% było retroelementami Gypsy (Fig. 1B i Supplemental Table S19). Datowanie molekularne oszacowało, że ta szczególna klasa retrotranspozonów LTR uległa amplifikacji ~38 Mya, co jest zgodne z wypiętrzeniem Płaskowyżu Qinghai-Tybetańskiego (Fig. 1C).

Szybka ewolucja rodzin genów związanych z patogennością grzybów

Jedną z najbardziej uderzających cech genomu O. sinensis jest brak wysoce homologicznych par genów. Spośród przewidywanych 7,939 genów kodujących białka, żadna para nie dzieliła >90% identyczności aminokwasów w sekwencjach kodujących i była tylko jedna para, która dzieliła >80% identyczności aminokwasów (Figura 2A i Supplemental Table S21). Cechę tę zaobserwowano również u blisko spokrewnionego C. militaris i ektomikoryzowego grzyba Tuber melanosporum 23. W porównaniu z innymi entomopatogenicznymi grzybami, takimi jak B. bassiana i C. militaris, rodziny wielogenowe u O. sinensis były nieliczne i stanowiły tylko 8,7% przewidywanego proteomu; większość rodzin genów miała tylko dwóch członków (Supplemental Figure S12). Tempo zyskiwania genów było uderzająco niższe niż utraty genów, a wśród 7 800 rodzin genów znalezionych w najnowszym wspólnym przodku (MRCA) Hypocreales, 1 756 zostało pozornie utraconych w O. sinensis (Fig. 2B). Tak zwarta przestrzeń kodowania genów genomu O. sinensis sugeruje naturę tego wysoce wyspecjalizowanego grzyba o niskiej zdolności adaptacji do wielu wskazówek środowiskowych.

Aby zrozumieć ewolucję rodzin genów, które odnoszą się do patogeniczności grzybów i adaptacji do surowych środowisk, zbadaliśmy właściwości funkcjonalne rodzin genów, które uległy ekspansji lub skurczeniu w O. sinensis. Genom O. sinensis wykazał znaczną ekspansję rodzin genów, które są głównie zaangażowane w patogenność grzybów, w tym aktywność peroksydazy (PF01328; P < 0.01), hydrolazę serynową (PF03959; P < 0,01), deuterolizynową metaloproteazę (M35) peptydazę (PF02102; P < 0,01) i cytochrom P450 (PF00067; P < 0,01) (Supplemental Table S23). Co ciekawe, odkryliśmy, że rozszerzone rodziny genów są również funkcjonalnie wzbogacone w kategorię Pfam oksydoreduktazy glukoza-metanol-cholina (GMC) zaangażowanej w metabolizm ecdysteroidów w moltingu u owadów (Supplemental Table S23). W porównaniach z innymi grzybami entomopatogennymi, ekspansja rodziny genów w linii O. sinensis była również obserwowana z nadreprezentacją terminów Pfam związanych z adaptacją do niskiej temperatury (PF06772; P < 0,01) (Supplemental Table S23).

Aby uniknąć infekcji patogenów grzybowych, gospodarze owadów często szybko wytwarzają mnóstwo reaktywnych form tlenu (ROS), aby bezpośrednio zabić patogeny. W odpowiedzi patogeny wykształciły podczas ewolucji system obrony antyoksydacyjnej ROS, którego peroksydazy, działające jako enzymy wymiatające ROS, są uważane za jeden z najbardziej znaczących i integralnych składników24, 25. Wśród ekspandowanych genów u O. sinensis aktywność peroksydaz była jedną z wysoko wzbogaconych kategorii funkcjonalnych (Supplemental Table S23). Przeszukiwanie za pomocą ukrytego modelu Markowa (HMM) ujawniło 42 (~0,53%) geny peroksydazy u O. sinensis, których liczba była znacznie większa niż u C. militaris (28) i drożdży (21), sugerując, że dwukrotna ekspansja genów peroksydazy może potencjalnie skutkować silną zdolnością do wspomagania detoksykacji ROS u O. sinensis (Fig. 3A i Supplemental Table S26). Wśród tych 42 genów peroksydazy, haloperoksydaza (haem) jest najbardziej obfita, stanowiąc ~16,67% wszystkich wykrytych peroksydaz (Rys. 3B). W przeciwieństwie do innych blisko spokrewnionych gatunków grzybów, które całkowicie pozbawione są typowej peroksydoksyny 2-Cysteinowej, O. sinensis nadal zachowuje jedną kopię (Rys. 3B). Wykazano wcześniej, że peroksydoksyna 2-Cys odgrywa rolę w odpowiedzi na różne poziomy stresu oksydacyjnego u Vibrio vulnificus 26. Analiza porównawcza wykazała, że zachowany gen u O. sinensis należy do typu Prx1, o którym donoszono, że jest funkcjonalnie konserwowany27 i ulega ekspresji tylko wtedy, gdy komórki są narażone na wysokie poziomy H2O2 generowane egzogennie26.

W przeciwieństwie do mechanizmu infekcji patogenów roślinnych (PPs), który wymaga enzymów węglowodanowo-czynnych (CAZymes) do degradacji złożonej ściany komórkowej rośliny28, patogeny owadzie (IPs) zazwyczaj infekują swoich gospodarzy poprzez penetrację kutykuli29. Aby to sprawdzić, porównaliśmy O. sinensis i cztery kolejne patogeny owadzie (M. anisopliae, M. acridum, C. militaris i B. bassiana) z czterema patogenami roślin (Fusarium graminearum, Magnaporthe grisea, Grosmannia clavigera i Botrytis cinerea) (Supplemental Table S17). Nasze wyniki wykazały, że patogeny owadów miały więcej proteaz (średnio 362 w IPs vs. 342 w PPs; P < 0,43) i kinaz białkowych (średnio 151 w IPs vs. 119 w PPs; P < 0,0014) do degradacji kutykuli owadów w porównaniu z patogenami roślinnymi (Tabele uzupełniające S27-29). W przeciwieństwie do tego, patogeny roślinne zawierały więcej CAZymes niż patogeny owadów do degradacji ściany komórkowej roślin (średnio 161 w IPs vs. 231 w PPs) (Tabele uzupełniające S30-32). Wyłączając patogeny roślinne, O. sinensis miał mniej genów kodujących proteazy (260) niż inne patogeny owadzie, takie jak M. anisopliae (437), M. acridum (361), C. militaris (355) i B. bassiana (396). Jednakże ~35% tych proteaz scharakteryzowanych u O. sinensis zawiera peptyd sygnałowy, który najprawdopodobniej był zaangażowany w interakcje patogen-gospodarz (tabele uzupełniające S10 i S34), co jest wartością większą niż u innych grzybów entomopatogennych (średnio 20%). Podobnie jak w przypadku innych patogenów owadzich, kilka rodzin celulaz, w tym GH7, GH45 i GH51, również zmniejszyło się lub było nieobecnych u O. sinensis (Tabela S30).

Badaliśmy również profile ekspresji genów w trzech stadiach rozwojowych O. sinensis, gdzie stosunek długości grzyba do owada osiągał ~1,20×, ~1,75× i ~2,20×. Wyniki pokazują, że w sumie 411 genów ulegało różnej ekspresji (DEG) pomiędzy trzema stadiami rozwojowymi (Supplemental Figure S14). Funkcjonalna anotacja tych 411 DEG wykazała, że były one głównie zaangażowane w patogenność grzybów, takie jak rodzina 16 hydrolaz glikozylowych (PF00722; FDR < 0.01), cytochrom P450 (PF00067; FDR < 0.01) i superrodzina głównych mediatorów (PF07690; FDR < 0.05). Ponadto, funkcjonalnie wzbogacone były również geny kodujące enzymy związane z mitochondrialnym łańcuchem oddechowym, takie jak rodzina NAD-zależnych epimerazy/dehydrataz (PF01370; FDR < 0,01) i N-końcowa domena BCS1 (PF08740; FDR < 0,01) (Supplementary Table S33).

Pozytywna darwinowska selekcja służy jako siły napędowe dla patogenności grzybów

Pozytywna selekcja niewątpliwie odegrała krytyczną rolę w ewolucji różnorodnych organizmów żyjących w środowiskach wysokogórskich Płaskowyżu Qinghai-Tybetańskiego, a wiele cech fenotypowych prawdopodobnie wykazuje takie sygnatury selekcji3,4,5. Spośród 1499 ortologów o wysokiej wiarygodności pojedynczej kopii, wspólnych dla O. sinensis i innych 12 gatunków grzybów, 163 pozytywnie wyselekcjonowane geny (PSGs) zostały zidentyfikowane w O. sinensis za pomocą testu współczynnika prawdopodobieństwa (LRT; P < 0,05) (Supplemental Table S35). Spośród nich, jeden gen (OSIN3929; tutaj nazwany OsPRX1) został funkcjonalnie zaangażowany w aktywność peroksydazy (Fig. 3C). Gen ten jest członkiem rodziny peroksiredoksyn z 1-cysteiną i jest wysoce homologiczny do PRX1 (YBL064C) w S. cerevisiae 30. Geny PRX1 w S. cerevisiae i dwóch patogennych dla człowieka grzybach, A. fumigatus i C. albicans, są funkcjonalnie konserwowane i wymagane do detoksykacji wybuchu oksydacyjnego w komórkach gospodarza31, 32. W szczególności, delecja PRX1 w znanym patogenie ryżu, Magnaporthe oryzae, spowodowała prawie całkowitą utratę patogenności, co sugeruje, że ta peroksydaza jest kluczowa dla interakcji gospodarz-patogen27. Co znamienne, kilka genów zaangażowanych w interakcje gospodarz-patogen, w tym biogeneza peroksysomalna, kinaza białkowa i metalopeptydazy, również zostały wykryte jako podlegające pozytywnej selekcji (Ryc. 3C).

Ewolucja systemu kojarzenia

W grzybach ascomycetous, system kojarzenia jest zwykle kontrolowany przez locus typu kojarzenia (MAT)33. Nasza analiza sekwencjonowania genomu wykazała, że O. sinensis nie tylko posiadał gen typu kojarzenia MAT1-2-1 w obrębie idiomorfu MAT1-2, ale także miał trzy geny typu kojarzenia (tj. MAT1-1-1, MAT1-1-2 i MAT1-1-3) w obrębie idiomorfu MAT1-1 (Supplemental Figure S15B). Cecha ta została zweryfikowana przy użyciu resekwencjonowania całego genomu 31 naturalnych populacji z niemal całego zasięgu geograficznego, wskazując, że O. sinensis jest rzeczywiście homotaliczny (Supplemental Figure S15A i Supplemental Table S36). Cecha ta jest skrajnie różna od jej blisko spokrewnionych patogenów grzybowych, takich jak Tolypocladium inflatum (MAT1-2)22, C. militaris (MAT1-1)20, B. bassiana (MAT1-1)21, M. anisopliae (MAT1-1)19 i M. acridum (MAT1-2)19, które są heterotaliczne i posiadają tylko jeden locus typu kojarzenia. Podobnie jak w przypadku dobrze znanego homotalicznego patogenu roślinnego Fusarium graminearum 34, organizacja tych dwóch loci MAT u O. sinensis ujawniła status fuzji w obrębie idiomorficznego regionu genomowego, który był szczególnie wzbogacony w retrotranspony LTR. Rozszczepienie między homotalicznym O. sinensis i heterotalicznym C. militaris oszacowano na blisko 174,2 MYA (Supplemental Figure S13C) i podlegał on wielokrotnym przekształceniom systemu kojarzenia z samozgodnego na samozgodny w trakcie swojej historii ewolucyjnej (Supplemental Figure S15C), przypominając nitkowaty ascomycete z rodzaju Neurospora 35.

Dywergencja populacji w oparciu o szerokości geograficzne na Płaskowyżu Qinghai-Tybetańskim

Aby zbadać relacje w całym genomie i dywergencję populacji, zebraliśmy i resekwencjonowaliśmy 31 dostępów O. sinensis w całym jego znanym zakresie dystrybucji, w tym prowincje Qinghai, Sichuan, Yunnan i Gansu oraz Tybetański Region Autonomiczny na Płaskowyżu Qinghai-Tybetańskim (Supplemental Figure S16 i Supplemental Table S37). Wygenerowaliśmy w sumie 183 miliony sparowanych odczytów (~36,68 Gb sekwencji) o średniej głębokości ~10,1× (dane surowe) (Supplemental Table S38). Z tych danych wygenerowaliśmy zestaw 816 960 polimorfizmów pojedynczego nukleotydu (SNPs) i 48 092 ścisłych indeli (insercji i delecji) w celu oceny pokrewieństwa między populacjami O. sinensis (Supplemental Figures S18-19 i Supplemental Table S39). Większość wariantów genomowych (71,1%) została zmapowana do regionów intergenicznych z podzbiorem mapującym do regionów kodujących (23,3% składającym się z 101 997 synonimicznych i 88 224 nonsynonimicznych SNP o współczynniku substytucji 0,86) (Supplemental Figure S19 i Supplemental Table S39). Drzewo filogenetyczne skonstruowane przy użyciu zestawów danych SNP podzieliło 31 dostępów na trzy geograficznie odrębne grupy od regionów położonych na niskich i wysokich szerokościach geograficznych (Ryc. 4A), co zostało wzmocnione przez PCA przy użyciu pierwszego i drugiego wektora własnego (Ryc. 4B i Supplemental Figure S20A). Zmienna liczba domniemanych populacji przodków (K) wykazała, że gdy K = 3, trzy odrębne grupy są zgodne z PCA i rekonstrukcją filogenetyczną (Fig. 4C i Supplemental Figure S20B). Niektóre akcesje z grupy nisko położonej wykazują silne dowody domieszki i są bardziej rozproszone w porównaniu z pozostałymi dwiema grupami, wskazując na większą różnorodność genetyczną, prawdopodobnie ze względu na wspólne polimorfizmy przodków i/lub niedawne zdarzenia introgresji (Fig. 4D,E). Oszacowana statystyka zróżnicowania populacji (F ST ) wśród tych trzech grup opartych na szerokości geograficznej dodatkowo ujawniła bazalny charakter regionu nizinnego, w szczególności populacji z dystryktu Nyingchi w Tybecie, co zostało dodatkowo potwierdzone przez znacznie podwyższone zróżnicowanie nukleotydów (π) w obrębie grupy i obniżone zróżnicowanie populacji z pozostałymi dwiema grupami wysokogórskimi. (Rys. 4C-F).

Dalej badaliśmy geny, na które wpływ miały różne poziomy zawartości SNP i mutacje niesynonimiczne (Supplemental Table S40). Analiza wzbogacenia funkcjonalnego 100 genów o najwyższej zawartości SNP i/lub mutacji niesynonimicznych wykazała, że są one głównie zaangażowane w metabolizm grzybowych metabolitów wtórnych, takich jak dehydrataza syntazy polyketydów (PF14765; FDR < 0,01), domena KR (PF08659; FDR < 0,01) i domena kondensacji (PF00668; FDR < 0,01). Wzbogacone zostały również kategorie funkcjonalne związane z biosyntezą kwasów tłuszczowych, takie jak domena transferazy acylowej (PF00698; FDR < 0,01) i syntazy beta-ketoacylu (PF00109 i PF02801; FDR < 0,01) (tabele uzupełniające S41-42).