Genom szekvenálás, összeszerelés és annotáció

Az O. sinensis-t a kínai Tibet Nyingchi körzetéből szekvenáltuk. WGS-elemzést végeztünk a Roche 454 és az Illumina HiSeq 2000 következő generációs szekvenáló platformokkal. Ez ~5,4 Gb nagyságú tiszta szekvencia-adatsorokat eredményezett, így körülbelül 45,1-szeres genomlefedettséget kaptunk (Kiegészítő S1 táblázat). A genom méretét ~124,08 Mb-ra és ~119,8 Mb-ra becsültük az áramlási citometria és a szekvenált leolvasások 17-mer mélységeloszlása alapján (Kiegészítő S1-2. ábra és Kiegészítő S2. táblázat). Az O. sinensis genomját először Roche 454 hosszú olvasatokból állítottuk össze a Newbler16 segítségével, majd az előzetesen összeállított kontigokat Illumina párpáros szekvenáló olvasatokkal állványoztuk össze az SSPACE17 segítségével. Ez végül egy ~116,4 Mb méretű genom-összeállítást eredményezett, amely a becsült genom méretének ~97%-át fedi le, és 156 scaffoldot (>2 Kb) tartalmaz ~3 Mb ScafN50 értékkel és 9 141 kontigot (N50 = 21 423 bp) (1. táblázat és kiegészítő táblázatok S3-4). A genom-összeállítás minőségének validálása érdekében először összehangoltuk az O. sinensis összes, a nyilvános adatbázisokban elérhető DNS-ét és kifejezett szekvenciacímkéjét (EST), és 98,85%-os, illetve 95,33%-os leképezési arányt kaptunk (Supplemental Table S5). Másodszor, az összes tiszta Roche 454 hosszú olvasatot (~1,84 Gb) leképeztük az összeállított genomszekvenciákra, ami 99,01%-os leképezési arány mellett szinte tökéletes illeszkedést mutatott (Supplemental Table S5). Harmadszor, az általunk összeállított transzkriptek jó illeszkedést mutattak az összeállított genomhoz; a 11 742 transzkript 91,29%-át sikerült leképezni (transzkript-lefedettség ≥80% és azonosság ≥90%; Supplemental Table S5). Végül a BUSCO18 segítségével értékeltük az O. sinensis összeszerelésünk teljességét; az 1 315 várható Ascomycota BUSCO konzervált gén 94,0%-át, 4,0%-át és 1,8%-át teljesnek, töredezettnek és hiányzónak azonosítottuk az O. sinensis összeszerelésünkben (Supplemental Table S5).

A génjóslás segítésére ~15,05 Gb RNS-szekvenálási (RNS-Seq) adatot generáltunk, amelyet összesen hat, a három fő fejlődési stádiumot reprezentáló könyvtárból nyertünk (S4. kiegészítő ábra és S6,7. kiegészítő táblázatok). Ab initio predikcióval, fehérje- és EST-illesztésekkel, EvidenceModeler fésüléssel és további szűréssel kombinálva 7939 fehérjekódoló gént határoztunk meg (1. táblázat és Supplemental Table S8). Ezen előrejelzett gének megközelítőleg 97,0%-át, illetve 71,51%-át sikerült funkcionálisan besorolni és RNS-Seq adatokkal alátámasztani (Kiegészítő táblázatok S9-11). Az Ascomycota vonalból származó BUSCO-t használva megállapítottuk továbbá, hogy a gének 94,4%-a teljes, 3,6%-a, 1,8%-a és 0,2%-a teljes, töredékes, hiányzó és duplikált volt, ami a génannotációnk jó minőségét jelzi (Supplemental Table S11). Homológkeresést is végeztünk, és annotáltuk a nem kódoló RNS (ncRNS) géneket, így 146 transzfer RNS (tRNS) gént, 33 riboszomális RNS (rRNS) gént, 70 kis nukleoláris RNS (snoRNS) gént és 15 kis nukleáris RNS (snRNS) gént találtunk (Supplemental Figure S6 és Supplemental Table S12). Az ismétlődő szekvenciák annotációja bemutatta, hogy a transzponálható elemek (TE-k) az összeszerelt genom körülbelül 74,67%-át és a nyers olvasatok 80,07%-át tették ki, ami azt jelzi, hogy az össze nem rakott genom ~5,45%-a TE-kből áll (Kiegészítő táblázatok S13-14). A GC-tartalom 43,09% volt a teljes genomban és 61,49% a kódoló szekvenciákban (Kiegészítő S3 ábra; Kiegészítő S4 és S8 táblázatok). Annotáltunk 8918 egyszerű szekvencia ismétlődést, amelyek értékes genetikai markereket biztosítanak a kínai hernyógomba jövőbeli nemesítési programjainak támogatására (Kiegészítő táblázatok S15-16 és Kiegészítő ábra S7).

Retrotranszpozonok által vezérelt genom-expanzió és a nem kollineáris gének tömeges eltávolítása

A genomméretek összehasonlítása azt mutatta, hogy az O. sinensis genomja közel 3,4-szer nagyobb, mint a Hypocreales családba tartozó más entomopatogén gombáké (Supplemental Table S17 és Supplemental Figure S8A). Az ismétlődő szekvenciaelemzés kimutatta, hogy ez a terjeszkedés elsősorban a transzpozíciós elemek gyors elszaporodásának köszönhető. Az O. sinensis genom-összeállításának mintegy 74,67%-a ismétlődő szekvenciákból állt (Kiegészítő táblázatok S13-14), ami kivételesen nagyobb, mint a Metarhizium anisopliae (~0,98%)19, Metarhizium acridum (~1,52%)19, Cordyceps militaris (~3,04%)20 és Beauveria bassiana (~2,03%)21 esetében (P < 4,822e-07) (Kiegészítő ábra S8B). A MULE elemek voltak a leggyakoribbak, az O. sinensis genom ~1,6%-át (~1,9 Mb) és az O. sinensis DNS-transzpozonok több mint 59%-át tették ki. A retrotranszpozonok, főként a hosszú terminális ismétlődésű (LTR) retrotranszpozonok az O. sinensis genom ~59,76%-át tették ki, és nagymértékű elszaporodásuk körülbelül ~38 millió évvel ezelőtt (MYA) következett be (Supplemental Figure S9).

Az LTR retrotranszpozonok gyors amplifikációjával ellentétben, amely az O. sinensis genom terjeszkedését hajtotta, egy másik figyelemre méltó jellemző a fehérjekódoló gének drámai elvesztése az O. sinensis vonalban más entomopatogén gombákhoz képest. Az O. sinensis összesen 7939 fehérjekódoló génjéhez képest más entomopatogén gombákban, pl. Metarhizium anisopliae (10 582)19, Metarhizium acridum (9 849)19, Cordyceps militaris (9 684)20, Beauveria bassiana (10 366)21 és Tolypocladium inflatum (9 998)22 átlagosan több mint 10 095 gén volt (1. táblázat). A génszám ilyen mértékű csökkenését a nem kollineáris gének azonosítása és az O. sinensis és a C. militaris genomjai közötti szinteniablokkok összehasonlító elemzése is bizonyította. Összesen 308 szintenikus blokkot azonosítottunk, amelyek a C. militaris genom közel 72,7%-át (~23,4 Mb a C. militarisban vs. ~43,5 Mb az O. sinensisben) fedik le (1A ábra; S10 kiegészítő ábra és S18-19 kiegészítő táblázatok). E szintenikus genomi régiók közül a nem kollineáris gének száma csökkent az O. sinensisben (2127) a C. militarishoz képest (3259), de az ismétlődő szekvenciák száma nőtt (23,8 Mb az O. sinensisben vs. 0,40 Mb a C. militarisban) (1B. ábra és S19. kiegészítő táblázat). Az O. sinensisben elveszett 2468 gén funkcionális annotációja azt mutatta, hogy ezek főként az aminosav-anyagcserében vesznek részt, például az aminosavak bioszintézisében (ko01230), az arginin- és prolin-anyagcserében (ko00330) és a tirozin-anyagcserében (ko00350) (S11. kiegészítő ábra és S20. kiegészítő táblázat). Figyelemre méltó, hogy az ismétlődő szekvenciák közel 81%-a ebben a 308 szintenikus blokkban LTR retrotranszpozonok voltak, amelyek 40,4%-a Gypsy retroelem (1B ábra és Supplemental Table S19). A molekuláris kormeghatározás becslése szerint az LTR retrotranszpozonoknak ez a bizonyos osztálya ~38 Mya alatt amplifikálódott, ami összhangban van a Qinghai-Tibeti-fennsík kiemelkedésével (1C ábra).

A gombák patogenitásával kapcsolatos géncsaládok gyors evolúciója

Az O. sinensis genom egyik legszembetűnőbb jellemzője a magasan homológ génpárok hiánya. Az előre jelzett 7939 fehérjekódoló gén közül egyetlen pár sem rendelkezett >90%-os aminosavazonossággal a kódoló szekvenciákban, és csak egy olyan pár volt, amely >80%-os aminosavazonossággal rendelkezett (2A ábra és Supplemental Table S21). Ez a jellegzetesség a közeli rokon C. militaris és az ektomikorrhiza gomba Tuber melanosporum 23 esetében is megfigyelhető volt. Más entomopatogén gombákkal, például a B. bassianával és a C. militarisszal összehasonlítva az O. sinensisben a multigéncsaládok száma korlátozott volt, és a prediktált proteomnak csak 8,7%-át tették ki; a legtöbb géncsaládnak csak két tagja volt (S12. kiegészítő ábra). A géngyarapodás aránya feltűnően alacsonyabb volt, mint a génvesztésé, és a Hypocreales legfrissebb közös ősében (MRCA) talált 7800 géncsalád közül 1756 látszólag elveszett az O. sinensisben (2B ábra). Az O. sinensis genomjának ilyen kompakt génkódolási tere arra utal, hogy ez az erősen specializált gomba kevéssé képes alkalmazkodni a többféle környezeti jelenséghez.

A gombák patogenitásával és a zord környezethez való magashegyi alkalmazkodással kapcsolatos géncsaládok evolúciójának megértése érdekében megvizsgáltuk az O. sinensisben bővülésen vagy szűkülésen átesett géncsaládok funkcionális tulajdonságait. Az O. sinensis genomban jelentős bővülést mutattak a géncsaládok, amelyek főként a gombák patogenitásában vesznek részt, beleértve a peroxidáz aktivitást (PF01328; P < 0.01), szerin-hidroláz (PF03959; P < 0,01), deuterolysin metalloproteáz (M35) peptidáz (PF02102; P < 0,01) és citokróm P450 (PF00067; P < 0,01) (S23 kiegészítő táblázat). Érdekes módon azt találtuk, hogy a kibővített géncsaládok funkcionálisan feldúsultak a rovarok vedlésénél az ecdiszteroid-anyagcserében részt vevő glükóz-metanol-kolin (GMC) oxidoreduktáz Pfam-kategóriában is (Supplemental Table S23). Más entomopatogén gombákkal való összehasonlításban az O. sinensis vonalban is megfigyelhető volt a géncsaládok bővülése az alacsony hőmérséklethez való alkalmazkodással feltételezhetően összefüggő Pfam-terminusok felülreprezentáltságával (PF06772; P < 0,01) (Supplemental Table S23).

Ezzel szemben az összehúzódási állapotot mutató géncsaládok elsősorban a transzportfolyamatokban és az energiaanyagcserében vettek részt, mint például az ABC-transzporterek (PF00005; P < 0,01), az aminosavpermeáz (PF00324; P < 0,01) és az ATP-szintáz (PF00306; P < 0,05) (Supplemental Table S24). E géncsaládok dinamikus evolúciója mellett további 1077 (~13,57%) fajspecifikus gént detektáltunk az O. sinensisben (2C ábra). Ezek közül 318 (~29,53%) gén funkcionálisan annotálható volt, és szignifikánsan feldúsultak a keményítőkötéssel (GO: 2001070; P < 0,01), patogenezissel (GO: 0009405; P < 0.01), és a sejtfal makromolekula katabolikus folyamat (GO: 0016998; P < 0,01) (Supplemental Table S25).

A gombás kórokozók fertőzésének elkerülése érdekében a rovargazdák gyakran gyorsan rengeteg reaktív oxigén species-t (ROS) termelnek a kórokozók közvetlen elpusztítása érdekében. Válaszként a kórokozók az evolúció során kifejlesztették a ROS-antioxidáns védelmi rendszert, amelynek egyik kiemelkedő és szerves komponensének a ROS-mentő enzimként működő peroxidázokat tekintik24, 25 . Az O. sinensisben kibővített gének között a peroxidáz-aktivitás volt az egyik erősen feldúsult funkcionális kategória (Supplemental Table S23). A rejtett Markov-modell (HMM) keresések 42 (~0,53%) peroxidáz gént mutattak ki az O. sinensisben, amelyek száma feltűnően nagyobb volt, mint a C. militarisban (28) és az élesztőben (21), ami arra utal, hogy a peroxidáz gének kétszeres bővülése potenciálisan erős ROS-mérgezést segítő képességet eredményezhet az O. sinensisben (3A ábra és Supplemental Table S26). E 42 peroxidáz gén közül a haloperoxidáz (haem) a leggyakoribb, amely az összes kimutatott peroxidáz ~16,67%-át teszi ki (3B ábra). Más, közeli rokon gombafajokkal ellentétben, amelyekből teljesen hiányzik a tipikus 2-cisztein peroxiredoxin, az O. sinensis még mindig megtartott egy példányt (3B ábra). Korábban kimutatták, hogy a 2-Cys-peroxiredoxin szerepet játszik a Vibrio vulnificusban a különböző szintű oxidatív stresszre adott válaszreakciókban 26. Egy összehasonlító elemzés kimutatta, hogy az O. sinensisben megtartott gén a Prx1 típusba tartozik, amelyről azt jelentették, hogy funkcionálisan konzervált27 és csak akkor fejeződik ki, ha a sejtek exogén módon előállított magas H2O2-szintnek vannak kitéve26.

A növényi kórokozók (PP-k) fertőzési mechanizmusával ellentétben, amely a komplex növényi sejtfal lebontásához szénhidrát-aktív enzimeket (CAZimek) igényel28, a rovarkórokozók (IP-k) jellemzően a kutikula áthatolásával fertőzik gazdaszervezetüket29. Ennek tesztelésére összehasonlítottuk az O. sinensis és négy további rovarkórokozó (M. anisopliae, M. acridum, C. militaris és B. bassiana) a négy növényi kórokozóval (Fusarium graminearum, Magnaporthe grisea, Grosmannia clavigera és Botrytis cinerea) (S17. kiegészítő táblázat). Eredményeink azt mutatták, hogy a rovarkórokozók több proteázzal (átlagosan 362 az IP-kben vs. 342 a PP-kben; P < 0,43) és proteinkinázzal (átlagosan 151 az IP-kben vs. 119 a PP-kben; P < 0,0014) bontották le a rovarkutikulát a növényi kórokozókhoz képest (Kiegészítő táblázatok S27-29). Ezzel szemben a növényi kórokozók a növényi sejtfal lebontásához több CAZymet tartalmaztak, mint a rovarkórokozók (átlagosan 161 az IP-kben vs. 231 a PP-kben) (Kiegészítő táblázatok S30-32). A növényi kórokozókat leszámítva az O. sinensis figyelemre méltó módon kevesebb proteázt kódoló génnel rendelkezett (260), mint más rovarkórokozók, például a M. anisopliae (437), M. acridum (361), C. militaris (355) és B. bassiana (396). Azonban ezeknek az O. sinensisben jellemzett proteázoknak ~35%-a tartalmaz olyan szignálpeptidet, amely nagyobb valószínűséggel vett részt a patogén-gazda interakciókban (Kiegészítő táblázatok S10 és S34), ami nagyobb, mint más entomopatogén gombákban (átlagosan 20%). A többi rovarkórokozóhoz hasonlóan az O. sinensisben is több cellulázcsalád, köztük a GH7, GH45 és GH51, csökkent vagy hiányzott (Kiegészítő táblázat S30).

Az O. sinensis három fejlődési stádiumában is vizsgáltuk a génexpressziós profilokat, a gomba vs. rovar hosszarányok ~1,20×, ~1,75× és ~2,20× elérésével. Az eredmények azt mutatják, hogy összesen 411 gén volt differenciálisan expresszált (DEG) a három fejlődési stádium között (Supplemental Figure S14). E 411 DEG funkcionális annotációja azt mutatta, hogy főként a gombák patogenitásában vesznek részt, mint például a 16-os glikozil-hidroláz család (PF00722; FDR < 0,01), a citokróm P450 (PF00067; FDR < 0,01) és a major facilitátor szupercsalád (PF07690; FDR < 0,05). Emellett a mitokondriális légzési lánccal kapcsolatos enzimeket kódoló gének is funkcionálisan gazdagodtak, mint például a NAD-függő epimeráz/dehidrátáz család (PF01370; FDR < 0,01) és a BCS1 N-terminális domén (PF08740; FDR < 0,01) (S33. kiegészítő táblázat).

A pozitív darwini szelekció a gombák patogenitásának hajtóerejeként szolgál

A pozitív szelekció kétségtelenül döntő szerepet játszott a Qinghai-Tibeti-fennsík magaslati környezetben élő változatos szervezetek evolúciójában, és számos fenotípusos tulajdonság valószínűleg ilyen szelekciós jegyeket mutat3,4,5. Az O. sinensis és a többi 12 gombafaj között közös, nagy megbízhatóságú egykópiás ortológok közül 163 pozitívan szelektált gént (PSG) azonosítottunk az O. sinensisben az elágazási hely valószínűségi arányteszt (LRT; P < 0,05) segítségével (Supplemental Table S35). Közülük egy gén (OSIN3929; itt OsPRX1 néven) funkcionálisan a peroxidáz-aktivitásban játszott szerepet (3C ábra). Ez a gén az 1-cisztein peroxiredoxin család tagja, és nagymértékben homológ az S. cerevisiae 30 PRX1 (YBL064C) génjével. A PRX1 gének az S. cerevisiae-ben és két humán patogén gombában, az A. fumigatusban és a C. albicansban funkcionálisan konzerváltak és szükségesek a gazdasejteken belüli oxidatív robbanás detoxifikációjához31, 32 . Különösen a PRX1 deléciója a jól ismert rizspatogénben, a Magnaporthe oryzae-ben a patogenitás majdnem teljes elvesztését eredményezte, ami arra utal, hogy ez a peroxidáz kulcsfontosságú a gazdatest-patogén kölcsönhatásokban27. Feltűnő, hogy számos, a gazdapatogén kölcsönhatásokban szerepet játszó gén, köztük a peroxiszomális biogenezis, a protein kináz és a metallopeptidázok is pozitív szelekció alatt álltak (3C. ábra).

Párzási rendszer evolúciója

Aszkomycetes gombákban a párzási rendszert általában a párzási típus (MAT) lókusz irányítja33. A genomszekvenálási elemzésünk szerint az O. sinensis nemcsak a MAT1-2-2 idiomorfon belül rendelkezett a MAT1-2 párosodási típus génnel, hanem a MAT1-1 idiomorfon belül három párosodási típus génnel (azaz MAT1-1-1, MAT1-1-2 és MAT1-1-3) is (Supplemental Figure S15B). Ezt a jellemzőt 31 természetes populáció teljes genom-újraszekvenálásával igazolták szinte a teljes földrajzi elterjedési területről, ami azt jelzi, hogy az O. sinensis valóban homotallikus (S15A kiegészítő ábra és S36 kiegészítő táblázat). Ez a tulajdonság rendkívül különbözik a közeli rokon gombakórokozóktól, mint például a Tolypocladium inflatum (MAT1-2)22, C. militaris (MAT1-1)20, B. bassiana (MAT1-1)21, M. anisopliae (MAT1-1)19 és M. acridum (MAT1-2)19, amelyek heterotallikusak és csak egyetlen párosodási lókusszal rendelkeznek. A jól ismert homotallikus növényi kórokozóhoz, a Fusarium graminearumhoz34 hasonlóan az O. sinensisben e két MAT-lókusz szerveződése feltárta az idiomorf genomikus régión belüli fúziós státuszt, amely különösen gazdag volt LTR retrotranszpononokban. A homotallikus O. sinensis és a heterotallikus C. militaris közötti szétválást közel 174,2 MYA-ra becsülték (Kiegészítő ábra S13C), és az evolúciós történelem során a párosodási rendszer többszörös átalakulásának volt kitéve öninkompatibilisről önkompatibilisre (Kiegészítő ábra S15C), hasonlóan a Neurospora 35 filamentózus aszkomycéta nemzetséghez.

Populációs divergencia a Qinghai-Tibeti-fennsík szélességi fokai alapján

A genomszintű kapcsolatok és a populációs divergencia vizsgálatához az O. sinensis 31 akcesszióját gyűjtöttük össze és szekvenáltuk újra az ismert elterjedési területéről, beleértve Qinghai, Sichuan, Yunnan és Gansu tartományokat, valamint a Qinghai-Tibeti-fennsík Tibeti Autonóm Területét (S16. kiegészítő ábra és S37. kiegészítő táblázat). Összesen 183 millió párosított végű olvasatot generáltunk (~36,68 Gb szekvencia), ~10,1×-es átlagos mélységgel (nyers adatok) (Supplemental Table S38). Ezekből az adatokból 816 960 egynukleotid-polimorfizmust (SNP) és 48 092 szigorú indelt (beillesztéseket és deléciókat) generáltunk az O. sinensis populációk közötti rokonsági viszonyok felméréséhez (S18-19. kiegészítő ábra és S39. kiegészítő táblázat). A genomiális variánsok többsége (71,1%) az intergenikus régiókba került leképezésre, egy részhalmaz pedig a kódoló régiókba (23,3%, amely 101 997 szinonim és 88 224 nem szinonim SNP-ből állt, 0,86-os helyettesítési aránnyal) (S19. kiegészítő ábra és S39. kiegészítő táblázat). Az SNP-adatkészletek felhasználásával szerkesztett filogenetikai fa a 31 hozzáférést három, földrajzilag elkülönülő csoportba osztotta, az alacsony szélességi és a magas szélességi régiók között (4A. ábra) – ezt a megállapítást az első és második sajátvektort használó PCA is megerősítette (4B. ábra és S20A kiegészítő ábra). A feltételezett őspopulációk számának (K) változtatása azt mutatta, hogy ha K = 3, akkor a három különböző csoport összhangban van a PCA és a filogenetikai rekonstrukcióval (4C. ábra és S20B kiegészítő ábra). Az alacsony szélességű csoport néhány hozzáférése erős bizonyítékot mutat a keveredésre, és a másik két csoporthoz képest jobban szétszóródott, ami nagyobb genetikai diverzitásra utal, ami valószínűleg a közös ősi polimorfizmusoknak és/vagy a közelmúltbeli introgressziós eseményeknek köszönhető (4D,E ábra). A becsült populációs differenciálódási statisztika (F ST ) e három szélességi fokon alapuló csoport között tovább mutatta az alacsony szélességi fokú régió bazális jellegét, különösen a tibeti Nyingchi körzetből származó populációkét, amit a csoporton belüli jelentősen megnövekedett nukleotiddiverzitás (π) és a másik két magas szélességi fokú csoporttal való alacsonyabb populációs differenciálódás is alátámaszt. (4C-F ábra).

Megvizsgáltuk továbbá a különböző szintű SNP-tartalmak és nem szinonim mutációk által érintett géneket (Supplemental Table S40). A legmagasabb SNP-tartalommal és/vagy nem szinonim mutációkkal rendelkező 100 gén funkcionális gazdagodási elemzése azt mutatja, hogy ezek főként a gombák másodlagos metabolitjainak anyagcseréjében vesznek részt, mint például a poliketid-szintáz dehidratáz (PF14765; FDR < 0,01), KR-domén (PF08659; FDR < 0,01) és kondenzációs domén (PF00668; FDR < 0,01). A zsírsav bioszintézishez kapcsolódó funkcionális kategóriák is feldúsultak, mint például az aciltranszferáz domén (PF00698; FDR < 0,01) és a béta-ketoacil szintáz (PF00109 és PF02801; FDR < 0,01) (S41-42. kiegészítő táblázatok).