Genoomsequencing, -assemblage en -annotatie

We hebben O. sinensis gesequeneerd uit het Nyingchi District van Tibet, China. We voerden een WGS analyse uit met de next-generation sequencing Roche 454 en Illumina HiSeq 2000 platforms. Dit genereerde schone sequentie datasets van ~ 5,4 Gb, waardoor ongeveer 45,1-voudige genoom dekking, respectievelijk (Supplementele Tabel S1). We schatten dat het genoom grootte is ~ 124.08 Mb en ~ 119.8 Mb op basis van flowcytometrie en 17-mer diepteverdeling van de sequenced leest, respectievelijk (Supplementary Figures S1-2 en Supplemental Table S2). De O. sinensis genoom werd eerst geassembleerd uit Roche 454 lange leest met behulp van Newbler 16, gevolgd door scaffolding voorgemonteerde contigs met Illumina mate paar sequencing leest met behulp van SSPACE 17. Dit leverde uiteindelijk een ~ 116.4 Mb genoom-assemblage die ~ 97% van de geschatte genoom-grootte dekt, en bevat 156 scaffolds (> 2 Kb) met een ScafN50 waarde van ~ 3 Mb en 9.141 contigs (N50 = 21.423 bp) (tabel 1 en supplementaire tabellen S3-4). Om de kwaliteit van de genoom assemblage te valideren, hebben we eerst alle DNA en express sequentie tags (ESTs) van O. sinensis uitgelijnd die beschikbaar waren in de publieke databases en karteringspercentages van 98,85% en 95,33% verkregen, respectievelijk (Supplementele Tabel S5). Ten tweede, we gekoppeld alle schone Roche 454 lange leest (~ 1,84 Gb) aan het geassembleerde genoom sequenties, met een bijna perfecte uitlijning met een mapping rate van 99,01% (Supplementary Table S5). Ten derde, de transcripten die we geassembleerd toonde een goede uitlijning met het geassembleerde genoom; van 11.742 transcripten, 91,29% werden in kaart gebracht (transcript dekking ≥80% en identiteit ≥90%; Supplementary Table S5). Tenslotte evalueerden we de volledigheid van onze O. sinensis assemblage met behulp van BUSCO18; 94,0%, 4,0% en 1,8% van de 1.315 verwachte Ascomycota BUSCO geconserveerde genen werden geïdentificeerd als compleet, gefragmenteerd en ontbrekend in onze O. sinensis assemblage, respectievelijk (Supplementele Tabel S5).

We genereerden ~15,05 Gb van RNA-sequencing (RNA-Seq) gegevens verkregen uit een totaal van zes bibliotheken die de drie belangrijkste ontwikkelingsstadia vertegenwoordigen om genvoorspelling te helpen (Supplementaire figuur S4 en Supplementaire tabellen S6,7). In combinatie met ab initio voorspelling, eiwit en EST alignments, EvidenceModeler combineren en verder filteren, definieerden we 7.939 eiwit-coderende genen (tabel 1 en supplementaire tabel S8). Van deze voorspelde genen kon ongeveer 97,0% en 71,51% functioneel worden geclassificeerd en ondersteund door RNA-Seq gegevens, respectievelijk (supplementaire tabellen S9-11). Met behulp van de BUSCO van Ascomycota lineage, vonden we verder dat 94,4%, 3,6%, 1,8% en 0,2% van de genen respectievelijk compleet, gefragmenteerd, ontbrekend en gedupliceerd waren, wat duidt op een goede kwaliteit van onze genannotatie (Supplementele Tabel S11). We hebben ook homologe zoekacties uitgevoerd en niet-coderende RNA (ncRNA) genen geannoteerd, wat 146 transfer RNA (tRNA) genen, 33 ribosomaal RNA (rRNA) genen, 70 kleine nucleolaire RNA (snoRNA) genen, en 15 kleine nucleaire RNA (snRNA) genen opleverde (Supplementaire Figuur S6 en Supplementaire Tabel S12). De annotatie van herhaalde sequenties bleek dat transposable elementen (TEs) goed voor ongeveer 74,67% van het geassembleerde genoom en 80,07% van de ruwe leest, wat aangeeft dat ~ 5,45% van het niet-geassembleerde genoom bestaat uit TEs (Supplementaire tabellen S13-14). Het GC-gehalte was 43,09% in het hele genoom en 61,49% in de coderende sequenties (supplementair Figuur S3; supplementaire tabellen S4 en S8). We annoteerden 8.918 simple sequence repeats die waardevolle genetische markers zullen opleveren om toekomstige veredelingsprogramma’s van Chinese rupsschimmel te ondersteunen (Supplementaire Tabellen S15-16 en Supplementaire Figuur S7).

Retrotransposon-gedreven genoom-expansie en massale verwijdering van niet-collineaire genen

Vergelijking van de genoomgrootte toonde aan dat het genoom van O. sinensis bijna 3,4-maal groter was dan dat van andere entomopathogene schimmels in de Hypocreales-familie (Supplementele Tabel S17 en Supplementele Figuur S8A). Analyse van de herhaalde sequentie toonde aan dat deze uitbreiding vooral te wijten was aan een snelle proliferatie van transponeerbare elementen. Ongeveer 74.67% van het genoom van O. sinensis bestond uit herhalingssequenties (supplementaire tabellen S13-14), uitzonderlijk groter dan die gerapporteerd in Metarhizium anisopliae (~0.98%)19, Metarhizium acridum (~1.52%)19, Cordyceps militaris (~3.04%)20 en Beauveria bassiana (~2.03%)21 (P < 4.822e-07) (supplementaire figuur S8B). De MULE elementen waren met name het meest talrijk, goed voor ~1.6% (~1.9 Mb) van het O. sinensis genoom en meer dan 59% van de DNA transposons in O. sinensis. Retrotransposons, meestal long-terminal repeat (LTR) retrotransposons, omvatten ~59.76% van het O. sinensis genoom en grootschalige proliferatie van die ongeveer op ~38 miljoen jaar geleden (MYA) plaatsvond (Supplemental Figure S9).

In tegenstelling tot de snelle amplificatie van LTR retrotransposons die de expansie van het O. sinensis genoom aansturen, is een ander opmerkelijk kenmerk het dramatische verlies van eiwit-coderende genen in de O. sinensis lineage in vergelijking met andere entomopathogene schimmels. Vergeleken met een totaal van 7.939 proteïne-coderende genen in O. sinensis waren er gemiddeld meer dan 10.095 genen in andere entomopathogene schimmels, b.v. Metarhizium anisopliae (10.582)19, Metarhizium acridum (9.849)19, Cordyceps militaris (9.684)20, Beauveria bassiana (10.366)21 en Tolypocladium inflatum (9.998)22 (Tabel 1). Een dergelijke vermindering van het gennummer werd verder aangetoond door de identificatie van niet-collineaire genen en vergelijkende analyse van syntenieblokken tussen de genomen van O. sinensis en C. militaris. We identificeerden een totaal van 308 syntenische blokken die bijna 72.7% (~23.4 Mb in C. militaris vs. ~43.5 Mb in O. sinensis) van het C. militaris genoom omspannen (Fig. 1A; Supplemental Figure S10 en Supplemental Tables S18-19). Van deze syntenische genomische regio’s was er een afname van niet-collineaire genen in O. sinensis (2.127) vergeleken met C. militaris (3.259) maar een toename van herhaalde sequenties (23.8 Mb in O. sinensis vs. 0.40 Mb in C. militaris) (Fig. 1B en Supplementary Table S19). Functionele annotatie van de 2.468 genen die verloren gingen in O. sinensis toonde aan dat ze voornamelijk betrokken waren bij het aminozuurmetabolisme, zoals biosynthese van aminozuren (ko01230), arginine en proline metabolisme (ko00330), en tyrosine metabolisme (ko00350) (Supplemental Figuur S11 en Supplemental Tabel S20). Opmerkelijk is dat bijna 81% van de herhaalde sequenties in deze 308 syntenische blokken LTR retrotransposons waren, waarvan 40,4% Gypsy retro-elementen waren (Fig. 1B en Supplementele Tabel S19). Moleculaire datering schatte dat deze specifieke klasse van LTR retrotransposons ~38 Mya amplificeerden, wat consistent is met de opheffing van het Qinghai-Tibetaans Plateau (Fig. 1C).

Snelle evolutie van genfamilies gerelateerd aan schimmelpathogeniciteit

Een van de meest opvallende kenmerken van het genoom van O. sinensis is de schaarste aan zeer homologe genenparen. Van de voorspelde 7.939 eiwit-coderende genen, deelde geen enkel paar >90% aminozuur-identiteit in coderende sequenties en er was slechts één paar dat >80% aminozuur-identiteit deelde (Fig. 2A en Supplementele Tabel S21). Dit kenmerk werd ook waargenomen bij de nauw verwante C. militaris en de ectomycorrhizaschimmel Tuber melanosporum 23. Vergeleken met andere entomopathogene schimmels zoals B. bassiana en C. militaris, waren multigene families in O. sinensis beperkt in aantal en maakten slechts 8.7% uit van het voorspelde proteoom; de meeste genfamilies hadden slechts twee leden (Supplementaire Figuur S12). De snelheid van genwinst was opvallend lager dan die van genverlies, en van de 7.800 genfamilies gevonden in de meest recente gemeenschappelijke voorouder (MRCA) van Hypocreales, gingen er 1.756 schijnbaar verloren in O. sinensis (Fig. 2B). Een dergelijke compacte gencoderingsruimte van het genoom van O. sinensis suggereert de aard van deze zeer gespecialiseerde schimmel met een geringe capaciteit om zich aan meerdere omgevingsfactoren aan te passen.

Om de evolutie te begrijpen van genfamilies die verband houden met schimmelpathogeniciteit en aanpassing aan de harde omgeving van het hoogland, hebben wij de functionele eigenschappen onderzocht van genfamilies die bij O. sinensis uitbreidingen of inkrimpingen hebben ondergaan. Het genoom van O. sinensis vertoonde een aanzienlijke uitbreiding van genfamilies die voornamelijk betrokken zijn bij schimmelpathogeniciteit, waaronder peroxidase activiteit (PF01328; P < 0.01), serine hydrolase (PF03959; P < 0.01), deuterolysine metalloprotease (M35) peptidase (PF02102; P < 0.01), en cytochroom P450 (PF00067; P < 0.01) (Supplementele Tabel S23). Interessant genoeg vonden we dat de uitgebreide genfamilies ook functioneel verrijkt zijn in de Pfam categorie van glucose-methanol-choline (GMC) oxidoreductase betrokken bij het ecdysteroïd metabolisme van de vervelling bij insecten (Supplementele Tabel S23). In vergelijkingen met andere entomopathogene schimmels, werd de genfamilie uitbreiding in de O. sinensis lineage ook waargenomen met over-vertegenwoordiging van Pfam termen putatief gerelateerd aan de aanpassing aan lage temperatuur (PF06772; P < 0.01) (Supplementele Tabel S23).

In tegenstelling daarmee waren genfamilies die een contractiestatus vertoonden voornamelijk betrokken bij het transportproces en energiemetabolisme, zoals ABC-transporteurs (PF00005; P < 0,01), aminozuurpermease (PF00324; P < 0,01), en ATP-synthase (PF00306; P < 0,05) (Supplemententabel S24). Afgezien van de dynamische evolutie van deze genfamilies, ontdekten we verder 1.077 (~13,57%) soortspecifieke genen in O. sinensis (Fig. 2C). Van hen konden 318 (~29,53%) genen functioneel geannoteerd worden en waren significant verrijkt in GO categorieën geassocieerd met zetmeelbinding (GO: 2001070; P < 0,01), pathogenese (GO: 0009405; P < 0.01), en celwand macromolecuul katabool proces (GO: 0016998; P < 0,01) (Supplementary Table S25).

Om de infectie van schimmelpathogenen te vermijden, produceren insecten gastheren vaak snel veel reactieve zuurstofsoorten (ROS) om pathogenen direct te doden. Als reactie hierop hebben ziekteverwekkers tijdens de evolutie het ROS antioxidant afweersysteem ontwikkeld, waarvan peroxidases, die fungeren als ROS-vangende enzymen, worden beschouwd als een van de meest prominente en integrale componenten24, 25. Onder de uitgebreide genen in O. sinensis, was peroxidase activiteit één van de hoogst verrijkte functionele categorieën (Supplementele Tabel S23). Zoeken met een verborgen Markov model (HMM) onthulde 42 (~0.53%) peroxidase genen in O. sinensis, waarvan het aantal opmerkelijk groter was dan dat van C. militaris (28) en gist (21), wat suggereert dat een tweevoudige uitbreiding van peroxidase genen mogelijk zou kunnen resulteren in een sterke capaciteit om ROS detoxificatie in O. sinensis te ondersteunen (Fig. 3A en Supplementele Tabel S26). Onder deze 42 peroxidase genen is haloperoxidase (haem) het meest overvloedig, goed voor ~16,67% van het totaal aan ontdekte peroxidases (Fig. 3B). In tegenstelling tot andere nauw verwante schimmelsoorten die het typische 2-Cysteine peroxiredoxine volledig missen, behoudt O. sinensis nog steeds één exemplaar (Fig. 3B). Eerder werd aangetoond dat 2-Cys peroxiredoxine een rol speelt in de reactie op verschillende niveaus van oxidatieve stress in Vibrio vulnificus 26. Een vergelijkende analyse toonde aan dat het behouden 2-Cys peroxiredoxine een rol speelt in de reactie op verschillende niveaus van oxidatieve stress. Een vergelijkende analyse toonde aan dat het behouden gen in O. sinensis tot het Prx1 type behoort, waarvan werd gerapporteerd dat het functioneel geconserveerd is27 en alleen tot expressie komt wanneer de cellen worden blootgesteld aan hoge niveaus van exogeen gegenereerd H2O2.

In tegenstelling tot het infectiemechanisme van plantpathogenen (PPs), waarbij koolhydraat-actieve enzymen (CAZymes) nodig zijn om de complexe plantencelwand af te breken28 , infecteren insectenpathogenen (IPs) hun gastheren doorgaans door de cuticula te penetreren29. Om dit te testen, vergeleken we O. sinensis met vier andere insectenpathogenen (M. anisopliae, M. acridum, C. militaris en B. bassiana) met de vier plantenpathogenen (Fusarium graminearum, Magnaporthe grisea, Grosmannia clavigera en Botrytis cinerea) (Supplemententabel S17). Onze resultaten toonden aan dat insectenpathogenen meer proteasen (gemiddeld 362 in IP’s vs. 342 in PP’s; P < 0,43) en proteïnekinasen (gemiddeld 151 in IP’s vs. 119 in PP’s; P < 0,0014) hadden om de insectenschub af te breken in vergelijking met plantpathogenen (Supplemententabellen S27-29). Daarentegen bezaten plantpathogenen meer CAZymen dan insectpathogenen voor de afbraak van plantencelwanden (gemiddeld 161 in IP’s vs. 231 in PP’s) (Supplemententabellen S30-32). Exclusief plantpathogenen, had O. sinensis opmerkelijk minder genen die codeerden voor proteasen (260) dan andere insectenpathogenen, zoals M. anisopliae (437), M. acridum (361), C. militaris (355), en B. bassiana (396). Echter, ~35% van deze proteasen gekarakteriseerd in O. sinensis bevatten een signaalpeptide waarvan het waarschijnlijker is dat het betrokken is geweest bij pathogeen-gastheer interacties (supplementaire tabellen S10 en S34), hetgeen groter is dan bij andere entomopathogene schimmels (gemiddeld 20%). Vergelijkbaar met de andere insectenpathogenen, namen verschillende cellulasefamilies, waaronder GH7, GH45, en GH51, ook af of waren afwezig in O. sinensis (Supplementele Tabel S30).

We onderzochten ook genexpressieprofielen in de drie ontwikkelingsstadia van O. sinensis met lengteverhoudingen van schimmel versus insect van ~1.20×, ~1.75× en ~2.20×. De resultaten tonen aan dat een totaal van 411 genen differentieel tot expressie kwamen (DEG) tussen de drie ontwikkelingsstadia (Supplementaire Figuur S14). Functionele annotatie van deze 411 DEG’s vond dat ze vooral betrokken waren bij schimmelpathogeniciteit, zoals glycosylhydrolasen familie 16 (PF00722; FDR < 0,01), cytochroom P450 (PF00067; FDR < 0,01) en major facilitator superfamilie (PF07690; FDR < 0,05). Bovendien waren genen die coderen voor enzymen die geassocieerd zijn met de mitochondriale ademhalingsketen ook functioneel verrijkt, zoals de NAD afhankelijke epimerase/dehydratase familie (PF01370; FDR < 0.01) en BCS1 N-terminaal domein (PF08740; FDR < 0.01) (Supplementary Table S33).

Positieve Darwinistische selectie dient als drijvende kracht voor schimmelpathogeniciteit

Positieve selectie heeft ongetwijfeld een cruciale rol gespeeld in de evolutie van diverse organismen die leven in hooggelegen milieus van het Qinghai-Tibetaans Plateau, en veel van de fenotypische kenmerken vertonen waarschijnlijk dergelijke selectiehandtekeningen3,4,5. Van de 1.499 hoog-vertrouwde single-copy orthologues gedeeld tussen O. sinensis en de andere 12 schimmelsoorten, werden 163 positief geselecteerde genen (PSGs) geïdentificeerd in O. sinensis door gebruik te maken van de branch-site likelihood ratio test (LRT; P < 0.05) (Supplementele Tabel S35). Van hen is één gen (OSIN3929; hier OsPRX1 genoemd) functioneel betrokken bij de peroxidase activiteit (Fig. 3C). Dit gen is lid van de peroxiredoxine familie met 1-cysteïne en zeer homoloog met PRX1 (YBL064C) in S. cerevisiae 30. PRX1-genen in S. cerevisiae en twee humane pathogene schimmels, A. fumigatus en C. albicans, zijn functioneel geconserveerd en vereist voor de ontgifting van de oxidatieve uitbarsting in gastheercellen31, 32. In het bijzonder leidde de deletie van PRX1 in de bekende rijstpathogeen Magnaporthe oryzae tot een bijna volledig verlies van pathogeniciteit, wat suggereert dat dit peroxidase een sleutelrol speelt in de interacties tussen gastheer en pathogeen27. Opvallend was dat verschillende genen die betrokken zijn bij gastheer-pathogeen interacties, waaronder peroxisomale biogenese, proteïnekinase en metallopeptidasen, ook onder positieve selectie bleken te staan (Fig. 3C).

Evolutie van het paringssysteem

In ascomycetente schimmels, wordt het paringssysteem gewoonlijk gecontroleerd door de mating-type (MAT) locus33. Onze genoomsequentieanalyse toonde aan dat O. sinensis niet alleen het MAT1-2-1 paringstype gen bezat binnen de MAT1-2 idiomorf, maar ook drie paringstype genen had (d.w.z., MAT1-1-1, MAT1-1-2, en MAT1-1-3) binnen de MAT1-1 idiomorf (Supplemental Figure S15B). Deze eigenschap werd geverifieerd met behulp van whole-genome resequencing van 31 natuurlijke populaties over bijna het gehele geografische verspreidingsgebied, wat aangeeft dat O. sinensis inderdaad homothallic is (Supplemental Figure S15A en Supplemental Table S36). Deze eigenschap verschilt sterk van de nauw verwante schimmelpathogenen, zoals Tolypocladium inflatum (MAT1-2)22, C. militaris (MAT1-1)20, B. bassiana (MAT1-1)21, M. anisopliae (MAT1-1)19, en M. acridum (MAT1-2)19, die heterothallic zijn en slechts een enkele paringstype locus bezitten. Vergelijkbaar met de bekende homothallische plantenpathogeen Fusarium graminearum 34, onthulde de organisatie van deze twee MAT loci in O. sinensis de fusiestatus binnen de idiomorfe genomische regio, die bijzonder verrijkt was in LTR retrotransponsons. De splitsing tussen de homothallische O. sinensis en heterothallische C. militaris werd geschat op bijna 174,2 MYA (Supplementaire Figuur S13C), en was onderhevig aan meerdere conversies van het paringssysteem van zelfincompatibel naar zelfcompatibel in de loop van zijn evolutionaire geschiedenis (Supplementaire Figuur S15C), lijkend op het filamenteuze ascomycetengeslacht Neurospora 35.

Populatie divergentie gebaseerd op breedtegraden op het Qinghai-Tibetaans Plateau

Om de genoom-brede relaties en populatie divergentie te onderzoeken, verzamelden en resequencing we 31 accessies van O. sinensis over zijn bekende verspreidingsgebied, inclusief de Qinghai, Sichuan, Yunnan en Gansu provincies en de Tibet Autonomous Region op het Qinghai-Tibetaans Plateau (Supplementaire Figuur S16 en Supplementaire Tabel S37). We genereerden een totaal van 183 miljoen gepaarde-end leest (~ 36,68 Gb van sequenties) met een gemiddelde diepte van ~ 10,1 × (ruwe gegevens) (Supplementary Table S38). Uit deze gegevens hebben we een set van 816.960 single-nucleotide polymorfismen (SNPs) en 48.092 strikte indels (inserties en deleties) gegenereerd om de verwantschap tussen populaties van O. sinensis te beoordelen (Supplementaire Figuren S18-19 en Supplementaire Tabel S39). De meerderheid van de genomische varianten (71,1%) werden in kaart gebracht in intergenische regio’s, met een subset die in kaart werden gebracht in de coderende regio’s (23,3% bestaande uit 101.997 synonieme en 88.224 niet-synonieme SNPs met een substitutie ratio van 0,86) (Supplementaire Figuur S19 en Supplementaire Tabel S39). De fylogenetische boom die met behulp van de SNP-datasets werd geconstrueerd, verdeelde de 31 toegangssoorten in drie geografisch gescheiden groepen, variërend van laag- tot hooggelegen gebieden (Fig. 4A) – een bevinding die werd versterkt door PCA met behulp van de eerste en tweede eigenvectoren (Fig. 4B en Supplemental Figure S20A). Het variëren van het aantal veronderstelde voorouderpopulaties (K) toonde aan dat wanneer K = 3, de drie verschillende groepen consistent zijn met de PCA en de fylogenetische reconstructie (Fig. 4C en Supplemental Figuur S20B). Sommige entiteiten van de groep op lage hoogte vertonen sterke aanwijzingen van vermenging en zijn meer verspreid dan de andere twee groepen, wat wijst op een grotere genetische diversiteit, mogelijk als gevolg van de gedeelde voorouderlijke polymorfismen en/of recente introgressie (Fig. 4D,E). De geschatte populatie-differentiatie statistiek (F ST ) tussen deze drie op breedtegraad gebaseerde groepen onthulde verder de basale aard van het laaggelegen gebied, met name populaties uit het Nyingchi district van Tibet, wat verder werd aangetoond door zijn aanzienlijk verhoogde nucleotide diversiteit (π) binnen de groep en verlaagde populatie-differentiatie met de andere twee hooggelegen groepen. (Fig. 4C-F).

We onderzochten verder de genen die beïnvloed werden door verschillende niveaus van SNP-inhoud en niet-synonieme mutaties (Supplementele Tabel S40). Functionele verrijkingsanalyse van de top 100 genen met de hoogste SNP-inhoud en/of niet-synonieme mutaties toont aan dat ze voornamelijk betrokken zijn bij het metabolisme van fungale secundaire metabolieten, zoals polyketide synthase dehydratase (PF14765; FDR < 0.01), KR-domein (PF08659; FDR < 0.01) en condensatiedomein (PF00668; FDR < 0.01). De functionele categorieën geassocieerd met vetzuur biosynthese waren ook verrijkt, zoals acyl transferase domein (PF00698; FDR < 0,01) en beta-ketoacyl synthase (PF00109 en PF02801; FDR < 0,01) (Supplementary Tables S41-42).