Sequenziamento del genoma, assemblaggio e annotazione

Abbiamo sequenziato O. sinensis dal distretto Nyingchi del Tibet, Cina. Abbiamo eseguito un’analisi WGS con le piattaforme di sequenziamento di nuova generazione Roche 454 e Illumina HiSeq 2000. Questo ha generato set di dati di sequenza pulita di ~ 5.4 Gb, producendo così circa 45.1-fold copertura del genoma, rispettivamente (Tabella supplementare S1). Abbiamo stimato che la dimensione del genoma è ~ 124,08 Mb e ~ 119,8 Mb sulla base della citometria a flusso e 17-mer distribuzione della profondità delle letture sequenziate, rispettivamente (figure supplementari S1-2 e tabella supplementare S2). Il genoma O. sinensis è stato prima assemblato da Roche 454 lunghe letture utilizzando Newbler 16, seguita da scaffolding contigs pre-assemblati con Illumina coppia mate sequenziamento legge utilizzando SSPACE 17. Questo finalmente prodotto un ~ 116,4 Mb genoma assembly che copre ~ 97% della dimensione stimata genoma, e contiene 156 scaffold (> 2 Kb) con un valore ScafN50 di ~ 3 Mb e 9.141 contigs (N50 = 21.423 bp) (Tabella 1 e Supplemento Tabelle S3-4). Per convalidare la qualità dell’assemblaggio del genoma, abbiamo prima allineato tutto il DNA e i tag di sequenza espressi (EST) di O. sinensis disponibili nei database pubblici e abbiamo ottenuto tassi di mappatura del 98,85% e 95,33%, rispettivamente (Tabella supplementare S5). In secondo luogo, abbiamo mappato tutte le letture pulite Roche 454 lunghe (~ 1,84 Gb) alle sequenze del genoma assemblate, mostrando un allineamento quasi perfetto con un tasso di mappatura del 99,01% (Tabella supplementare S5). In terzo luogo, i trascritti che abbiamo assemblato hanno mostrato un buon allineamento al genoma assemblato; di 11.742 trascrizioni, il 91,29% sono stati mappati (copertura trascrizione ≥80% e identità ≥90%; Tabella supplementare S5). Infine, abbiamo valutato la completezza del nostro assemblaggio O. sinensis utilizzando BUSCO18; 94,0%, 4,0% e 1,8% dei 1.315 geni conservati Ascomycota BUSCO previsti sono stati identificati come completi, frammentati e mancanti nel nostro assemblaggio O. sinensis, rispettivamente (Tabella supplementare S5).

Abbiamo generato ~15,05 Gb di dati di sequenziamento dell’RNA (RNA-Seq) ottenuti da un totale di sei librerie che rappresentano le tre principali fasi di sviluppo per aiutare la predizione del gene (Figura supplementare S4 e Tabelle supplementari S6,7). In combinazione con la predizione ab initio, allineamenti di proteine e EST, pettinatura EvidenceModeler e ulteriore filtraggio, abbiamo definito 7.939 geni codificanti proteine (Tabella 1 e Tabella supplementare S8). Di questi geni predetti, circa il 97,0% e il 71,51% potrebbe essere classificato funzionalmente e supportato da dati RNA-Seq, rispettivamente (Tabelle supplementari S9-11). Utilizzando il BUSCO dal lignaggio Ascomycota, abbiamo inoltre scoperto che il 94,4%, 3,6%, 1,8% e 0,2% dei geni erano completi, frammentati, mancanti e duplicati, rispettivamente, indicando una buona qualità della nostra annotazione genica (Tabella supplementare S11). Abbiamo anche eseguito ricerche di omologhi e annotato geni di RNA non codificanti (ncRNA), ottenendo 146 geni di RNA di trasferimento (tRNA), 33 geni di RNA ribosomiale (rRNA), 70 geni di piccoli RNA nucleolari (snoRNA) e 15 geni di piccoli RNA nucleari (snRNA) (Figura supplementare S6 e Tabella supplementare S12). L’annotazione delle sequenze ripetute ha presentato che gli elementi trasponibili (TEs) rappresentavano circa il 74.67% del genoma assemblato e il 80.07% delle letture grezze, indicando che ~ 5.45% del genoma non assemblato è costituito da TEs (tabelle supplementari S13-14). Il contenuto di GC era del 43,09% in tutto il genoma e del 61,49% nelle sequenze codificanti (Figura supplementare S3; Tabelle supplementari S4 e S8). Abbiamo annotato 8.918 ripetizioni di sequenza semplice che fornirà preziosi marcatori genetici per assistere i futuri programmi di allevamento del fungo cingolato cinese (Tabelle supplementari S15-16 e Figura supplementare S7).

Espansione del genoma guidata da retrotransposoni e rimozione massiccia di geni non collinari

Il confronto delle dimensioni del genoma ha mostrato che il genoma di O. sinensis era quasi 3,4 volte più grande di altri funghi entomopatogeni della famiglia Hypocreales (Tabella supplementare S17 e Figura supplementare S8A). L’analisi della sequenza di ripetizione ha rivelato che questa espansione era principalmente dovuta a una rapida proliferazione di elementi trasponibili. Circa il 74,67% del genoma di O. sinensis era composto da sequenze ripetute (Tabelle supplementari S13-14), eccezionalmente più grandi di quelle riportate in Metarhizium anisopliae (~0,98%)19, Metarhizium acridum (~1,52%)19, Cordyceps militaris (~3,04%)20 e Beauveria bassiana (~2,03%)21 (P < 4,822e-07) (Figura supplementare S8B). Gli elementi MULE erano in particolare il più abbondante, che rappresentano ~ 1,6% (~ 1,9 Mb) del genoma O. sinensis e più del 59% dei trasposoni DNA in O. sinensis. Retrotransposoni, per lo più retrotrasposoni a ripetizione long-terminale (LTR), comprendeva ~ 59,76% del genoma di O. sinensis e proliferazione su larga scala di cui si è verificato circa a ~ 38 milioni di anni fa (MYA) (Figura supplementare S9).

In contrasto con la rapida amplificazione dei retrotrasposoni LTR che guidano l’espansione del genoma di O. sinensis, un’altra caratteristica notevole è la drammatica perdita di geni codificanti proteine nel lignaggio di O. sinensis rispetto ad altri funghi entomopatogeni. Rispetto a un totale di 7.939 geni codificanti le proteine in O. sinensis c’erano più di 10.095 geni in media in altri funghi entomopatogeni, ad esempio Metarhizium anisopliae (10.582)19, Metarhizium acridum (9.849)19, Cordyceps militaris (9.684)20, Beauveria bassiana (10.366)21 e Tolypocladium inflatum (9.998)22 (Tabella 1). Tale riduzione del numero di geni è stata ulteriormente evidenziata dall’identificazione di geni non collinari e dall’analisi comparativa dei blocchi di sintenia tra i genomi di O. sinensis e C. militaris. Abbiamo identificato un totale di 308 blocchi sintenici che abbracciano quasi il 72,7% (~ 23,4 Mb in C. militaris vs. ~ 43,5 Mb in O. sinensis) del genoma C. militaris (Fig. 1A; Figura supplementare S10 e tabelle supplementari S18-19). Di queste regioni genomiche sinteniche, c’è stata una diminuzione dei geni non-collineari in O. sinensis (2.127) rispetto a C. militaris (3.259) ma un aumento delle sequenze ripetute (23,8 Mb in O. sinensis vs. 0,40 Mb in C. militaris) (Fig. 1B e Tabella supplementare S19). L’annotazione funzionale dei 2.468 geni che sono stati persi in O. sinensis ha mostrato che erano principalmente coinvolti nel metabolismo degli aminoacidi, come la biosintesi degli aminoacidi (ko01230), arginina e prolina metabolismo (ko00330), e tirosina metabolismo (ko00350) (Figura supplementare S11 e Tabella supplementare S20). In particolare, quasi 81% delle sequenze ripetute in questi 308 blocchi sintetici erano retrotrasposoni LTR, 40.4% dei quali erano retroelementi Gypsy (Fig. 1B e Tabella supplementare S19). Datazione molecolare stimato che questa particolare classe di retrotrasposoni LTR amplificato ~ 38 Mya, che è coerente con il sollevamento del Qinghai-Tibetan Plateau (Fig. 1C).

Rapida evoluzione delle famiglie di geni legati alla patogenicità dei funghi

Una delle caratteristiche più sorprendenti del genoma di O. sinensis è la penuria di coppie di geni altamente omologhi. Dei 7.939 geni codificanti proteine previsti, nessuna coppia condivideva >90% di identità aminoacidica nelle sequenze codificanti e c’era solo una coppia che condivideva >80% di identità aminoacidica (Fig. 2A e Tabella supplementare S21). Questa caratteristica è stata osservata anche nel C. militaris strettamente correlato e nel fungo ectomycorrhizal Tuber melanosporum 23. Rispetto ad altri funghi entomopatogeni come B. bassiana e C. militaris, famiglie multigene in O. sinensis erano limitati nel numero e comprendeva solo 8.7% del proteoma previsto; maggior parte delle famiglie di geni aveva solo due membri (Figura supplementare S12). Il tasso di guadagno del gene era sorprendentemente inferiore a quello della perdita del gene, e tra le 7.800 famiglie di geni trovate nell’antenato comune più recente (MRCA) delle Hypocreales, 1.756 sono state apparentemente perse in O. sinensis (Fig. 2B). Uno spazio di codifica genica così compatto del genoma di O. sinensis suggerisce la natura di questo fungo altamente specializzato con una bassa capacità di adattamento a molteplici spunti ambientali.

Per comprendere l’evoluzione delle famiglie di geni che si riferiscono alla patogenicità fungina e all’adattamento dell’altopiano ad ambienti difficili, abbiamo studiato le proprietà funzionali delle famiglie di geni che hanno subito espansioni o contrazioni in O. sinensis. Il genoma di O. sinensis ha mostrato una notevole espansione delle famiglie di geni che sono principalmente coinvolti nella patogenicità fungina, tra cui l’attività perossidasi (PF01328; P < 0.01), serina idrolasi (PF03959; P < 0.01), deuterolisina metalloproteasi (M35) peptidasi (PF02102; P < 0.01), e citocromo P450 (PF00067; P < 0.01) (Tabella supplementare S23). È interessante notare che le famiglie di geni espanse sono anche funzionalmente arricchiti nella categoria Pfam di glucosio-metanolo-colina (GMC) ossidoreduttasi coinvolti nel metabolismo ecdysteroid di stampo negli insetti (Tabella supplementare S23). Nel confronto con altri funghi entomopatogeni, l’espansione della famiglia di geni nel lignaggio O. sinensis è stato osservato anche con sovrarappresentare i termini Pfam putativamente correlati all’adattamento della bassa temperatura (PF06772; P < 0,01) (Supplemental Table S23).

Al contrario, le famiglie di geni che esibiscono lo stato di contrazione erano principalmente coinvolte nel processo di trasporto e nel metabolismo energetico, come i trasportatori ABC (PF00005; P < 0.01), la permeasi degli aminoacidi (PF00324; P < 0.01), e l’ATP sintasi (PF00306; P < 0.05) (Tabella supplementare S24). Oltre all’evoluzione dinamica di queste famiglie di geni, abbiamo ulteriormente rilevato 1.077 (~13,57%) geni specie-specifici in O. sinensis (Fig. 2C). Di loro, 318 (~ 29,53%) geni potrebbero essere annotati funzionalmente e sono stati significativamente arricchiti in categorie GO associate al legame dell’amido (GO: 2001070; P < 0,01), patogenesi (GO: 0009405; P < 0.01), e processo catabolico delle macromolecole della parete cellulare (GO: 0016998; P < 0.01) (Tabella supplementare S25).

Per evitare l’infezione di patogeni fungini, gli ospiti insetti spesso producono rapidamente molte specie reattive dell’ossigeno (ROS) per uccidere direttamente i patogeni. Come risposta, i patogeni hanno sviluppato il sistema di difesa antiossidante ROS durante l’evoluzione, di cui le perossidasi, che agiscono come enzimi di ROS-scavenging, sono considerate come uno dei componenti più importanti e integrali24, 25. Tra i geni espansi in O. sinensis, l’attività di perossidasi era una delle categorie funzionali altamente arricchite (Tabella supplementare S23). Hidden Markov model (HMM) ricerche rivelato 42 (~ 0,53%) geni perossidasi in O. sinensis, il cui numero è stato notevolmente più grande di quello di C. militaris (28) e lievito (21), suggerendo che una doppia espansione di geni perossidasi potrebbe potenzialmente risultare in una forte capacità di aiutare la detossificazione ROS in O. sinensis (Fig. 3A e Tabella supplementare S26). Tra questi 42 geni perossidasi, haloperoxidase (haem) è il più abbondante, che rappresenta ~ 16,67% del totale perossidasi rilevato (Fig. 3B). In contrasto con altre specie fungine strettamente correlate che mancano completamente della tipica perossidossina 2-Cisteina, O. sinensis ne conserva ancora una copia (Fig. 3B). 2-Cys peroxiredoxin è stato precedentemente dimostrato di svolgere un ruolo nel rispondere a diversi livelli di stress ossidativo in Vibrio vulnificus 26. Un’analisi comparativa ha rivelato che il gene conservato in O. sinensis appartiene al tipo Prx1, che è stato segnalato per essere funzionalmente conservato27 ed espresso solo quando le cellule sono esposte ad alti livelli di H2O2 generato esogenamente26.

A differenza del meccanismo di infezione dei patogeni delle piante (PPs), che richiede enzimi carboidrati-attivi (CAZymes) per degradare la parete cellulare complessa pianta28, insetto patogeni (IPs) tipicamente infettare i loro ospiti penetrando la cuticola29. Per verificare questo, abbiamo confrontato O. sinensis e altri quattro insetti patogeni (M. anisopliae, M. acridum, C. militaris e B. bassiana) con i quattro patogeni vegetali (Fusarium graminearum, Magnaporthe grisea, Grosmannia clavigera e Botrytis cinerea) (Tabella supplementare S17). I nostri risultati hanno dimostrato che gli insetti patogeni avevano più proteasi (in media, 362 in IPs vs. 342 in PPs; P < 0.43) e protein chinasi (in media, 151 in IPs vs. 119 in PPs; P < 0.0014) per degradare la cuticola insetto rispetto ai patogeni delle piante (Supplemental Tables S27-29). Al contrario, i patogeni delle piante hanno ospitato più CAZymes rispetto ai patogeni degli insetti per la degradazione della parete cellulare delle piante (in media, 161 in IPs contro 231 in PPs) (Tabelle supplementari S30-32). Escludendo i patogeni delle piante, O. sinensis aveva notevolmente meno geni che codificano proteasi (260) rispetto ad altri insetti patogeni, come M. anisopliae (437), M. acridum (361), C. militaris (355), e B. bassiana (396). Tuttavia, ~ 35% di queste proteasi caratterizzate in O. sinensis contengono un peptide di segnale che è più probabile essere stato coinvolto nelle interazioni patogeno-ospite (tabelle supplementari S10 e S34), che è più grande di quello in altri funghi entomopatogeni (in media, 20%). Simile agli altri insetti patogeni, diverse famiglie di cellulasi, tra cui GH7, GH45, e GH51, anche diminuito o erano assenti in O. sinensis (Tabella supplementare S30).

Abbiamo anche esaminato i profili di espressione genica attraverso le tre fasi di sviluppo di O. sinensis con rapporti di lunghezza del fungo vs insetto raggiungendo ~ 1.20×, ~ 1.75× e ~ 2.20×. I risultati mostrano che un totale di 411 geni erano differenzialmente espressi (DEG) tra le tre fasi di sviluppo (Figura supplementare S14). Annotazione funzionale di questi 411 DEG trovato che erano principalmente coinvolti nella patogenicità fungina, come glicosil idrolasi famiglia 16 (PF00722; FDR < 0.01), citocromo P450 (PF00067; FDR < 0.01) e superfamiglia facilitatore maggiore (PF07690; FDR < 0.05). Inoltre, i geni che codificano gli enzimi associati con la catena respiratoria mitocondriale erano anche funzionalmente arricchito, come NAD dipendente epimerasi/deidratasi famiglia (PF01370; FDR < 0.01) e BCS1 N-terminal domain (PF08740; FDR < 0.01) (Tabella supplementare S33).

Selezione positiva darwiniana serve come forze motrici per patogenicità fungina

Selezione positiva ha indubbiamente giocato un ruolo critico nell’evoluzione di diversi organismi che vivono in ambienti di alta quota del Qinghai-Tibetan Plateau, e molti dei tratti fenotipici sono suscettibili di mostrare tale selezione firme3,4,5. Dei 1.499 ortologhi a copia singola ad alta confidenza condivisi tra O. sinensis e le altre 12 specie fungine, 163 geni selezionati positivamente (PSG) sono stati identificati in O. sinensis utilizzando il branch-site likelihood ratio test (LRT; P < 0,05) (Tabella supplementare S35). Di questi, un gene (OSIN3929; qui chiamato OsPRX1) è stato implicato funzionalmente nell’attività della perossidasi (Fig. 3C). Questo gene è un membro della famiglia delle perossidossine con 1-cisteina e altamente omologo a PRX1 (YBL064C) in S. cerevisiae 30. I geni PRX1 in S. cerevisiae e in due funghi patogeni umani, A. fumigatus e C. albicans, sono funzionalmente conservati e richiesti per la detossificazione del burst ossidativo nelle cellule ospiti31, 32. In particolare, la delezione di PRX1 nel noto patogeno del riso, Magnaporthe oryzae, ha portato a una perdita quasi completa della patogenicità, suggerendo che questa perossidasi è fondamentale per le interazioni ospite-patogeno27. Sorprendentemente, diversi geni coinvolti nelle interazioni ospite-patogeno, tra cui biogenesi perossisomica, protein chinasi e metallopeptidasi, sono stati rilevati anche per essere sotto selezione positiva (Fig. 3C).

Evoluzione del sistema di accoppiamento

Nei funghi ascomiceti, il sistema di accoppiamento è solitamente controllato dal locus mating-type (MAT)33. La nostra analisi di sequenziamento del genoma trovato che O. sinensis non solo possedeva il gene MAT1-2-1 tipo di accoppiamento all’interno del MAT1-2 idiomorfo ma aveva anche tre geni mating-type (cioè, MAT1-1-1, MAT1-1-2 e MAT1-1-3) all’interno del MAT1-1 idiomorfo (Figura supplementare S15B). Questa caratteristica è stata verificata utilizzando il risequenziamento dell’intero genoma di 31 popolazioni naturali in quasi tutto l’intervallo geografico, indicando che O. sinensis è effettivamente omotallico (Figura supplementare S15A e Tabella supplementare S36). La caratteristica è estremamente diversa dai suoi patogeni fungini strettamente correlati, come Tolypocladium inflatum (MAT1-2)22, C. militaris (MAT1-1)20, B. bassiana (MAT1-1)21, M. anisopliae (MAT1-1)19, e M. acridum (MAT1-2)19, che sono eterotalli e possiedono solo un singolo locus mating-type. Simile al noto patogeno vegetale omotallico Fusarium graminearum 34, l’organizzazione di questi due loci MAT in O. sinensis ha rivelato lo stato di fusione all’interno della regione genomica idiomorfa, che era particolarmente arricchita in retrotraspononi LTR. La scissione tra l’omotallico O. sinensis e eterotallico C. militaris è stato stimato per verificarsi quasi 174,2 MYA (Figura supplementare S13C), ed è stato sottoposto a conversioni multiple del sistema di accoppiamento da auto-incompatibile a auto-compatibile nel corso della sua storia evolutiva (Figura supplementare S15C), simile al genere ascomicete filamentoso di Neurospora 35.

Divergenza delle popolazioni in base alle latitudini sull’altopiano Qinghai-Tibetano

Per esaminare le relazioni a livello di genoma e la divergenza delle popolazioni, abbiamo raccolto e risequenziato 31 accessioni di O. sinensis in tutto il suo range di distribuzione noto, comprese le province Qinghai, Sichuan, Yunnan e Gansu e la regione autonoma del Tibet sull’altopiano Qinghai-Tibetano (Figura supplementare S16 e Tabella supplementare S37). Abbiamo generato un totale di 183 milioni di letture paired-end (~ 36.68 Gb di sequenze) con una profondità media di ~ 10.1× (dati grezzi) (Tabella supplementare S38). Da questi dati, abbiamo generato un set di 816.960 polimorfismi a singolo nucleotide (SNPs) e 48.092 indel rigorosi (inserzioni e delezioni) per valutare la parentela tra le popolazioni di O. sinensis (figure supplementari S18-19 e tabella supplementare S39). La maggior parte delle varianti genomiche (71,1%) sono stati mappati alle regioni intergeniche con un sottoinsieme di mappatura alle regioni codificanti (23,3% composto da 101.997 sinonimi e 88.224 SNP non sinonimi con la razione di sostituzione di 0,86) (Figura supplementare S19 e Tabella supplementare S39). L’albero filogenetico costruito utilizzando i set di dati SNP ha diviso le 31 accessioni in tre gruppi geograficamente separati che vanno dalle regioni a bassa latitudine a quelle ad alta latitudine (Fig. 4A) – una scoperta che è stata rafforzata da PCA utilizzando il primo e il secondo autovettore (Fig. 4B e Figura supplementare S20A). Variando il numero di presunte popolazioni ancestrali (K) ha mostrato che quando K = 3, i tre gruppi distinti sono coerenti con la PCA e la ricostruzione filogenetica (Fig. 4C e Figura supplementare S20B). Alcune accessioni dal gruppo di bassa latitudine mostrano una forte evidenza della commistione e sono più dispersi rispetto agli altri due gruppi, indicando una maggiore diversità genetica forse a causa dei polimorfismi ancestrali condivisi e/o eventi di introgressione recenti (Fig. 4D,E). La stima della statistica di differenziazione della popolazione (F ST ) tra questi tre gruppi basati sulla latitudine ha rivelato ulteriormente la natura basale della regione di bassa latitudine, in particolare le popolazioni del distretto di Nyingchi del Tibet, che è stata ulteriormente evidenziata dalla sua diversità nucleotidica sostanzialmente elevata (π) all’interno del gruppo e la differenziazione della popolazione più bassa con gli altri due gruppi di alta latitudine. (Fig. 4C-F).

Abbiamo ulteriormente studiato i geni interessati da diversi livelli di contenuti SNP e mutazioni non sinonimi (Tabella supplementare S40). L’analisi di arricchimento funzionale dei primi 100 geni con il più alto contenuto di SNP e/o mutazioni non sinonime mostra che sono principalmente coinvolti nel metabolismo dei metaboliti secondari fungini, come la deidratasi polichetoide sintasi (PF14765; FDR < 0.01), il dominio KR (PF08659; FDR < 0.01) e il dominio di condensazione (PF00668; FDR < 0.01). Le categorie funzionali associate con la biosintesi degli acidi grassi sono stati anche arricchiti, come il dominio transferasi acile (PF00698; FDR < 0,01) e beta-chetoacil sintasi (PF00109 e PF02801; FDR < 0,01) (tabelle supplementari S41-42).