Sequenciação, montagem e anotação do genoma

Sequenciamos O. sinensis do Distrito Nyingchi do Tibete, China. Realizamos uma análise WGS com a sequenciação da próxima geração das plataformas Roche 454 e Illumina HiSeq 2000. Isto gerou conjuntos de dados de sequência limpa de ~5,4 Gb, produzindo assim aproximadamente 45,1 vezes a cobertura do genoma, respectivamente (Tabela Suplementar S1). Estimamos que o tamanho do genoma é de ~124.08 Mb e ~119.8 Mb com base na citometria de fluxo e distribuição de 17 metros de profundidade de leituras sequenciadas, respectivamente (Figuras Suplementares S1-2 e Tabela Suplementar S2). O genoma O. sinensis foi primeiramente montado a partir de leituras longas de Roche 454 usando Newbler16, seguido por contíguos pré-montados de andaimes com seqüenciamento de pares de Illumina mate com leituras usando SSPACE17. Isto finalmente produziu uma montagem do genoma de ~116.4 Mb que cobre ~97% do tamanho estimado do genoma, e contém 156 andaimes (>2 Kb) com um valor ScafN50 de ~3 Mb e 9.141 contíguos (N50 = 21.423 bp) (Tabela 1 e Tabelas Suplementares S3-4). Para validar a qualidade da montagem do genoma, primeiramente alinhamos todos os tags de DNA e sequências expressas (ESTs) de O. sinensis disponíveis nas bases de dados públicas e obtivemos taxas de mapeamento de 98,85% e 95,33%, respectivamente (Tabela Suplementar S5). Em segundo lugar, mapeamos todas as leituras limpas do Roche 454 (~1,84 Gb) para as seqüências do genoma montado, mostrando um alinhamento quase perfeito com uma taxa de mapeamento de 99,01% (Tabela Suplementar S5). Em terceiro lugar, as transcrições que montamos mostraram um bom alinhamento com o genoma montado; de 11.742 transcrições, 91,29% foram mapeadas (cobertura das transcrições ≥80% e identidade ≥90%; Tabela Suplementar S5). Finalmente, avaliamos a integridade da nossa montagem de O. sinensis usando BUSCO18; 94,0%, 4,0% e 1,8% dos 1.315 genes conservados de Ascomycota BUSCO esperados foram identificados como completos, fragmentados e ausentes em nossa montagem de O. sinensis, respectivamente (Tabela Suplementar S5).

Geramos ~15,05 Gb de dados de sequenciamento do RNA (RNA-Seq) obtidos de um total de seis bibliotecas representando os três principais estágios de desenvolvimento para auxiliar na predição do gene (Figura Suplementar S4 e Tabelas Suplementares S6,7). Em combinação com a previsão ab initio, os alinhamentos protéicos e EST, a penteação do EvidenceModeler e a filtragem posterior, definimos 7.939 genes codificadores de proteínas (Tabela 1 e Tabela Suplementar S8). Destes genes previstos, aproximadamente 97,0% e 71,51% poderiam ser classificados funcionalmente e suportados por dados do RNA-Seq, respectivamente (Tabelas Suplementares S9-11). Usando o BUSCO da linhagem Ascomycota, encontramos ainda que 94,4%, 3,6%, 1,8% e 0,2% dos genes estavam completos, fragmentados, ausentes e duplicados, respectivamente, indicando uma boa qualidade da nossa anotação gênica (Tabela Suplementar S11). Também realizamos pesquisas homólogas e anotamos genes não codificadores de RNA (ncRNA), produzindo 146 genes RNA de transferência (tRNA), 33 genes RNA ribossômico (rRNA), 70 pequenos genes RNA nucleolares (snoRNAs) e 15 pequenos genes RNA nucleares (snRNA) (Figura Suplementar S6 e Tabela Suplementar S12). A anotação de seqüências repetidas apresentou que elementos transponíveis (ETs) representaram aproximadamente 74,67% do genoma montado e 80,07% das leituras brutas, indicando que ~5,45% do genoma não montado consiste de ETs (Tabelas Suplementares S13-14). O conteúdo de GC foi 43,09% em todo o genoma e 61,49% em seqüências de codificação (Figura Suplementar S3; Tabelas Suplementares S4 e S8). Anotamos 8.918 repetições de sequências simples que fornecerão valiosos marcadores genéticos para auxiliar futuros programas de reprodução do fungo lagarta chinês (Tabelas Suplementares S15-16 e Figura Suplementar S7).

Expansão do genoma por retrotransposição e remoção maciça de genes não-colineares

Comparação do tamanho do genoma mostrou que o genoma O. sinensis era quase 3,4 vezes maior que outros fungos entomopatogênicos na família Hypocreales (Tabela Suplementar S17 e Figura Suplementar S8A). A análise da sequência repetida revelou que esta expansão se deveu principalmente a uma rápida proliferação de elementos transponíveis. Aproximadamente 74,67% do genoma do O. sinensis foi composto de sequências repetidas (Tabelas Suplementares S13-14), excepcionalmente maiores que as relatadas em 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 Suplementar S8B). Os elementos MULE foram notavelmente os mais abundantes, representando ~1,6% (~1,9 Mb) do genoma O. sinensis e mais de 59% dos transpositores de DNA em O. sinensis. Retrotransposões, na sua maioria retrotransposões de repetição a longo prazo (LTR), compreenderam ~59,76% do genoma O. sinensis e cuja proliferação em larga escala ocorreu aproximadamente a ~38 milhões de anos atrás (MYA) (Figura Suplementar S9).

Em contraste com a rápida amplificação dos retrotransposons LTR que impulsionaram a expansão do genoma O. sinensis, outra característica notável é a dramática perda de genes codificadores de proteínas na linhagem O. sinensis em comparação com outros fungos entomopatogênicos. Comparado a um total de 7.939 genes codificadores de proteína no O. sinensis, havia mais de 10.095 genes em média em outros fungos entomopatogênicos, por exemplo, 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 (Tabela 1). Tal redução do número de genes foi ainda evidenciada pela identificação de genes não-colineares e análise comparativa de bloqueios sintéticos entre os genomas O. sinensis e C. militaris. Nós identificamos um total de 308 blocos sintéticos que abrangem quase 72,7% (~23,4 Mb em C. militaris vs. ~43,5 Mb em O. sinensis) do genoma C. militaris (Fig. 1A; Figura Suplementar S10 e Tabelas Suplementares S18-19). Destas regiões do genoma sintético, houve uma diminuição de genes não-colineares em O. sinensis (2.127) comparado com C. militaris (3.259) mas um aumento de seqüências repetidas (23.8 Mb em O. sinensis vs. 0.40 Mb em C. militaris) (Fig. 1B e Tabela Suplementar S19). A anotação funcional dos 2.468 genes perdidos em O. sinensis mostrou que eles estavam principalmente envolvidos no metabolismo de aminoácidos, como biossíntese de aminoácidos (ko01230), metabolismo de arginina e prolina (ko00330), e metabolismo de tirosina (ko00350) (Figura Suplementar S11 e Tabela Suplementar S20). Notavelmente, quase 81% das sequências repetidas nesses 308 blocos sintéticos foram retrotransposições LTR, dos quais 40,4% foram retroelementos ciganos (Figura 1B e Tabela Suplementar S19). A datação molecular estimou que esta classe particular de retrotransposões LTR amplificou ~38 Mya, o que é consistente com a elevação do Planalto Qinghai-Tibetano (Fig. 1C).

Evolução rápida das famílias de genes relacionados à patogenicidade fúngica

Uma das características mais marcantes do genoma O. sinensis é a escassez de pares de genes altamente homólogos. Dos 7.939 genes codificadores de proteínas previstos, nenhum par compartilhou >90% de identidade de aminoácidos nas sequências de codificação e houve apenas um par que compartilhou >80% de identidade de aminoácidos (Fig. 2A e Tabela Suplementar S21). Esta característica também foi observada no C. militaris e no fungo ectomicorrízico Tuber melanosporum 23. Comparado com outros fungos entomopatogênicos como B. bassiana e C. militaris, as famílias multigenas em O. sinensis eram limitadas em número e compreendiam apenas 8,7% do proteoma previsto; a maioria das famílias de genes tinha apenas dois membros (Figura Suplementar S12). A taxa de ganho de genes foi surpreendentemente menor do que a de perda de genes, e entre as 7.800 famílias de genes encontradas no mais recente ancestral comum (MRCA) de Hypocreales, 1.756 foram aparentemente perdidas em O. sinensis (Fig. 2B). Um espaço de codificação gênica tão compacto do genoma O. sinensis sugere a natureza deste fungo altamente especializado com baixa capacidade de adaptação a múltiplos estímulos ambientais.

Para entender a evolução das famílias de genes que se relacionam com a patogenicidade fúngica e adaptação a ambientes agressivos, investigamos as propriedades funcionais de famílias de genes que sofreram expansões ou contrações em O. sinensis. O genoma O. sinensis mostrou uma expansão considerável das famílias de genes que estão principalmente envolvidas na patogenicidade fúngica, incluindo a atividade da peroxidase (PF01328; P < 0.01), hidrolase serina (PF03959; P < 0,01), deuterolysin metaloprotease (M35) peptidase (PF02102; P < 0,01), e citocromo P450 (PF00067; P < 0,01) (Tabela suplementar S23). Curiosamente, descobrimos que as famílias de genes expandidos também são enriquecidas funcionalmente na categoria Pfam de glucose-metanol-colina (GMC) oxidoreductase envolvida no metabolismo ecdisteróide de molting em insetos (Tabela Suplementar S23). Em comparações com outros fungos entomopatogênicos, a expansão da família genética na linhagem O. sinensis também foi observada com termos Pfam super-representativos putativamente relacionados à adaptação de baixa temperatura (PF06772; P < 0,01) (Tabela Suplementar S23).

Em contraste, famílias de genes exibindo estado de contração estiveram envolvidas principalmente no processo de transporte e metabolismo energético, tais como transportadores ABC (PF00005; P < 0,01), permease de aminoácidos (PF00324; P < 0,01), e ATP sintase (PF00306; P < 0,05) (Tabela Suplementar S24). Além da evolução dinâmica destas famílias de genes, detectamos ainda 1.077 (~13,57%) genes específicos de espécies em O. sinensis (Fig. 2C). Destes, 318 (~29,53%) genes puderam ser anotados funcionalmente e foram significativamente enriquecidos nas categorias GO associados à ligação do amido (GO: 2001070; P < 0,01), patogênese (GO: 0009405; P < 0.01), e processo catabólico de macromoléculas da parede celular (GO: 0016998; P < 0.01) (Tabela suplementar S25).

Para evitar a infecção de patógenos fúngicos, hospedeiros de insetos freqüentemente produzem rapidamente muitas espécies reativas de oxigênio (ROS) para matar diretamente os patógenos. Como resposta, os patógenos desenvolveram o sistema de defesa antioxidante ROS durante a evolução, dos quais as peroxidases, atuando como enzimas necrófagas ROS, são consideradas como um dos componentes mais proeminentes e integrais24, 25. Entre os genes expandidos em O. sinensis, a atividade da peroxidase foi uma das categorias funcionais altamente enriquecidas (Tabela Suplementar S23). Pesquisas do modelo Markov oculto (HMM) revelaram 42 (~0,53%) genes peroxidase no O. sinensis, cujo número foi notavelmente maior que o de C. militaris (28) e levedura (21), sugerindo que uma expansão dupla dos genes peroxidase pode potencialmente resultar em uma forte capacidade de auxiliar a desintoxicação ROS no O. sinensis (Fig. 3A e Tabela Suplementar S26). Entre estes 42 genes de peroxidase, a haloperoxidase (hemoglobina) é a mais abundante, representando ~16,67% do total de peroxidases detectadas (Fig. 3B). Ao contrário de outras espécies fúngicas intimamente relacionadas que carecem completamente da típica peroxiredoxina 2-Cysteine, O. sinensis ainda retém uma cópia (Fig. 3B). A peroxiredoxina 2-Cys demonstrou anteriormente ter um papel na resposta a diferentes níveis de stress oxidativo no Vibrio vulnificus 26. Uma análise comparativa revelou que o gene retido em O. sinensis pertence ao tipo Prx1, que foi reportado como funcionalmente conservado27 e expresso apenas quando as células são expostas a níveis elevados de H2O2 gerados exogenamente26.

Não se parece com o mecanismo de infecção dos patógenos vegetais (PPs), que requer enzimas carboidratos-ativas (CAZymes) para degradar a complexa parede celular das plantas28, patógenos de insetos (IPs) tipicamente infectam seus hospedeiros penetrando na cutícula29. Para testar isso, comparamos O. sinensis e mais quatro patógenos de insetos (M. anisopliae, M. acridum, C. militaris e B. bassiana) com os quatro patógenos vegetais (Fusarium graminearum, Magnaporthe grisea, Grosmannia clavigera e Botrytis cinerea) (Tabela suplementar S17). Nossos resultados demonstraram que os patógenos dos insetos tinham mais proteases (em média, 362 em IPs vs. 342 em PPs; P < 0,43) e proteínas kinases (em média, 151 em IPs vs. 119 em PPs; P < 0,0014) para degradar a cutícula do inseto em comparação com os patógenos vegetais (Tabelas Suplementares S27-29). Em contraste, os patógenos vegetais abrigavam mais CAZymes do que insetos patógenos para a degradação da parede celular das plantas (em média, 161 em IPs vs. 231 em PPs) (Tabelas Suplementares S30-32). Excluindo patógenos vegetais, O. sinensis tinha notavelmente menos genes codificando proteases (260) do que outros patógenos de insetos, tais como M. anisopliae (437), M. acridum (361), C. militaris (355), e B. bassiana (396). Entretanto, ~35% dessas proteases caracterizadas em O. sinensis contêm um peptídeo sinal que é mais provável que tenha estado envolvido em interações patogênico-hospitalares (Tabelas Suplementares S10 e S34), que é maior do que o de outros fungos entomopatogênicos (em média, 20%). Similar aos outros patógenos insetos, várias famílias de celulase, incluindo GH7, GH45 e GH51, também declinaram ou estiveram ausentes em O. sinensis (Tabela Suplementar S30).

Nós também examinamos perfis de expressão gênica através dos três estágios de desenvolvimento de O. sinensis com proporções de comprimento do fungo vs. insetos atingindo ~1,20×, ~1,75× e ~2,20×. Os resultados mostram que um total de 411 genes foram expressos diferentemente (DEG) entre os três estágios de desenvolvimento (Figura Suplementar S14). A anotação funcional desses 411 DEGs encontrou que eles estavam principalmente envolvidos na patogenicidade fúngica, como a família das glicosil hidrolases 16 (PF00722; FDR < 0,01), citocromo P450 (PF00067; FDR < 0,01) e superfamília facilitadora principal (PF07690; FDR < 0,05). Além disso, os genes que codificam as enzimas associadas à cadeia respiratória mitocondrial também foram enriquecidos funcionalmente, como a família de epimerase/deidratase dependente de NAD (PF01370; FDR < 0,01) e o domínio N-terminal BCS1 (PF08740; FDR < 0,01) (Tabela suplementar S33).

A seleção darwiniana positiva serve como força motriz para a patogenicidade fúngica

A seleção positiva tem sem dúvida desempenhado um papel crítico na evolução de diversos organismos que vivem em ambientes de alta altitude do Planalto Qinghai-Tibetano, e muitos dos traços fenotípicos são passíveis de mostrar tais assinaturas de seleção3,4,5. Dos 1.499 ortologues de alta confiança de cópia única compartilhados entre O. sinensis e as outras 12 espécies fúngicas, 163 genes selecionados positivamente (PSGs) foram identificados em O. sinensis usando o teste de razão de verossimilhança do ramo (LRT; P 0,05) (Tabela suplementar S35). Destes, um gene (OSIN3929; aqui denominado OsPRX1) foi implicado funcionalmente na atividade da peroxidase (Fig. 3C). Este gene é um membro da família da peroxiredoxina com 1-cisteína e altamente homólogo ao PRX1 (YBL064C) em S. cerevisiae 30. Os genes PRX1 em S. cerevisiae e dois fungos patogênicos humanos, A. fumigatus e C. albicans, são conservados funcionalmente e necessários para a desintoxicação da explosão oxidativa dentro das células hospedeiras31, 32. Em particular, a deleção de PRX1 no conhecido patógeno do arroz, Magnaporthe oryzae, resultou em uma perda quase completa da patogenicidade, sugerindo que esta peroxidase é a chave para as interações hospedeiro-patógeno27. Surpreendentemente, vários genes envolvidos nas interações hospedeiro-patógeno, incluindo a biogênese peroxisomal, proteína quinase e metallopeptidases, também foram detectados sob seleção positiva (Fig. 3C).

Evolução do sistema de acasalamento

Em fungos ascomicetos, o sistema de acasalamento é geralmente controlado pelo locus tipo acasalamento (MAT)33. Nossa análise de seqüenciamento genômico descobriu que O. sinensis não só possuía o gene do tipo acasalamento MAT1-2-1 dentro do idiomorfo MAT1-2, mas também tinha três genes do tipo acasalamento (isto é, MAT1-1-1, MAT1-1-2, e MAT1-1-3) dentro do idiomorfo MAT1-1 (Figura suplementar S15B). Esta característica foi verificada usando ressequenciamento de 31 populações naturais em quase toda a área geográfica, indicando que O. sinensis é de fato homoftálico (Figura Suplementar S15A e Tabela Suplementar S36). A característica é extremamente diferente de seus patógenos fúngicos estreitamente relacionados, tais como 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, que são heterotálicos e possuem apenas um único locus do tipo acasalamento. Similar ao conhecido patógeno vegetal homoftálico Fusarium graminearum 34, a organização destes dois MAT loci em O. sinensis revelou o estado de fusão dentro da região genômica idiomórfica, que foi particularmente enriquecida em retrotransponsos LTR. A divisão entre o O. sinensis homotálico e o C. militaris heterotálico foi estimada em quase 174,2 MYA (Figura Suplementar S13C), e foi submetido a múltiplas conversões do sistema de acasalamento de auto-incompatível para auto-compatível durante o curso de sua história evolutiva (Figura Suplementar S15C), assemelhando-se ao gênero ascomycete filamentoso de Neurospora 35.

Divergência populacional baseada nas latitudes do Planalto Qinghai-Tibetano

Para examinar as relações entre o genoma e a divergência populacional, coletamos e ressequenciamos 31 adesões de O. sinensis em toda a sua distribuição conhecida, incluindo as províncias de Qinghai, Sichuan, Yunnan e Gansu e a Região Autônoma Tibetana do Planalto Qinghai-Tibetano (Figura Suplementar S16 e Tabela Suplementar S37). Gerámos um total de 183 milhões de leituras em ponta de página (~36,68 Gb de sequências) com uma profundidade média de ~10,1× (dados brutos) (Tabela Suplementar S38). A partir desses dados, geramos um conjunto de 816.960 polimorfismos de nucleotídeos simples (SNPs) e 48.092 indels rigorosos (inserções e deleções) para avaliar a relação entre as populações de O. sinensis (Figuras Suplementares S18-19 e Tabela Suplementar S39). A maioria das variantes genômicas (71,1%) foi mapeada para regiões intergênicas com um subconjunto mapeado para as regiões codificadoras (23,3% consistindo de 101.997 sinônimos e 88.224 SNPs não-sinônimos com a ração de substituição de 0,86) (Figura Suplementar S19 e Tabela Suplementar S39). A árvore filogenética construída usando os conjuntos de dados do SNP dividiu as 31 adesões em três grupos geograficamente separados, variando de regiões de baixa latitude a regiões de alta latitude (Figura 4A) – um achado que foi reforçado pela APC usando o primeiro e segundo vetores próprios (Figura 4B e Figura Suplementar S20A). Variando o número de presumíveis populações ancestrais (K) mostrou que quando K = 3, os três grupos distintos são consistentes com a PCA e a reconstrução filogenética (Fig. 4C e Figura Suplementar S20B). Alguns acessos do grupo de baixa latitude apresentam forte evidência da mistura e estão mais dispersos em comparação com outros dois grupos, indicando maior diversidade genética possivelmente devido aos polimorfismos ancestrais compartilhados e/ou eventos de introgressão recentes (Fig. 4D,E). A estatística estimada de diferenciação populacional (F ST ) entre estes três grupos baseados na latitude revelou ainda a natureza basal da região de baixa latitude, particularmente as populações do distrito de Nyingchi no Tibete, o que foi ainda mais evidenciado pela sua diversidade de nucleotídeos substancialmente elevada (π) dentro do grupo e a menor diferenciação populacional com os outros dois grupos de alta latitude. (Fig. 4C-F).

Investigamos ainda mais os genes afetados por diferentes níveis de conteúdo de SNP e mutações não-sinônimas (Tabela suplementar S40). A análise do enriquecimento funcional dos 100 genes com maior conteúdo de SNP e/ou mutações não-sinônimas mostra que eles estão principalmente envolvidos no metabolismo de metabolitos secundários fúngicos, como a politase sintetase desidratase (PF14765; FDR < 0,01), domínio KR (PF08659; FDR < 0,01) e domínio de condensação (PF00668; FDR < 0,01). As categorias funcionais associadas à biossíntese de ácidos graxos também foram enriquecidas, tais como domínio acyl transferase (PF00698; FDR < 0,01) e beta-ketoacyl synthase (PF00109 e PF02801; FDR < 0,01) (Tabelas Suplementares S41-42).