Secuenciación, ensamblaje y anotación del genoma

Secuenciamos O. sinensis del distrito de Nyingchi del Tíbet, China. Realizamos un análisis WGS con las plataformas de secuenciación de próxima generación Roche 454 e Illumina HiSeq 2000. Esto generó conjuntos de datos de secuencias limpias de ~5,4 Gb, produciendo así una cobertura del genoma de aproximadamente 45,1 veces, respectivamente (Tabla Suplementaria S1). Estimamos que el tamaño del genoma es de ~124,08 Mb y ~119,8 Mb basándonos en la citometría de flujo y en la distribución de la profundidad de 17 marcadores de las lecturas secuenciadas, respectivamente (Figuras suplementarias S1-2 y Tabla suplementaria S2). El genoma de O. sinensis se ensambló primero a partir de las lecturas largas de Roche 454 utilizando Newbler16, seguido de un andamiaje de contigs preensamblados con las lecturas de secuenciación de pares de Illumina utilizando SSPACE17. Finalmente se obtuvo un ensamblaje del genoma de ~116,4 Mb que cubre el ~97% del tamaño estimado del genoma y contiene 156 scaffolds (>2 Kb) con un valor ScafN50 de ~3 Mb y 9.141 contigs (N50 = 21.423 bp) (Tabla 1 y Tablas Suplementarias S3-4). Para validar la calidad del ensamblaje del genoma, primero alineamos todo el ADN y las etiquetas de secuencia expresada (EST) de O. sinensis disponibles en las bases de datos públicas y obtuvimos tasas de mapeo del 98,85% y 95,33%, respectivamente (Tabla Suplementaria S5). En segundo lugar, mapeamos todas las lecturas largas limpias de Roche 454 (~1,84 Gb) con las secuencias genómicas ensambladas, mostrando una alineación casi perfecta con una tasa de mapeo del 99,01% (Tabla Suplementaria S5). En tercer lugar, los transcritos que ensamblamos mostraron una buena alineación con el genoma ensamblado; de 11.742 transcritos, el 91,29% fueron mapeados (cobertura del transcrito ≥80% e identidad ≥90%; Tabla Suplementaria S5). Por último, evaluamos la integridad de nuestro ensamblaje de O. sinensis utilizando BUSCO18; el 94,0%, el 4,0% y el 1,8% de los 1.315 genes conservados esperados de Ascomycota BUSCO se identificaron como completos, fragmentados y ausentes en nuestro ensamblaje de O. sinensis, respectivamente (Tabla Suplementaria S5).

Generamos ~15,05 Gb de datos de secuenciación de ARN (RNA-Seq) obtenidos a partir de un total de seis bibliotecas que representaban los tres principales estadios de desarrollo para ayudar a la predicción de genes (Figura Suplementaria S4 y Tablas Suplementarias S6,7). En combinación con la predicción ab initio, los alineamientos de proteínas y EST, el peinado EvidenceModeler y el filtrado posterior, definimos 7.939 genes codificadores de proteínas (Tabla 1 y Tabla Suplementaria S8). De estos genes predichos, aproximadamente el 97,0% y el 71,51% pudieron ser clasificados funcionalmente y apoyados por datos de RNA-Seq, respectivamente (Tablas Suplementarias S9-11). Utilizando el BUSCO del linaje Ascomycota, descubrimos además que el 94,4%, el 3,6%, el 1,8% y el 0,2% de los genes estaban completos, fragmentados, ausentes y duplicados, respectivamente, lo que indica una buena calidad de nuestra anotación genética (Tabla Suplementaria S11). También realizamos búsquedas de homólogos y anotamos genes de ARN no codificantes (ncRNA), obteniendo 146 genes de ARN de transferencia (tRNA), 33 genes de ARN ribosómico (rRNA), 70 genes de ARN nucleolar pequeño (snoRNAs) y 15 genes de ARN nuclear pequeño (snRNA) (Figura Suplementaria S6 y Tabla Suplementaria S12). La anotación de las secuencias repetidas presentó que los elementos transponibles (TEs) representaban aproximadamente el 74,67% del genoma ensamblado y el 80,07% de las lecturas crudas, indicando que ~5,45% del genoma no ensamblado consiste en TEs (Tablas Suplementarias S13-14). El contenido de GC fue del 43,09% en todo el genoma y del 61,49% en las secuencias codificantes (Figura Suplementaria S3; Tablas Suplementarias S4 y S8). Anotamos 8.918 repeticiones de secuencia simple que proporcionarán valiosos marcadores genéticos para ayudar a los futuros programas de mejora del hongo chino de la oruga (Tablas Suplementarias S15-16 y Figura Suplementaria S7).

Expansión del genoma impulsada por retrotransposones y eliminación masiva de genes no colineales

La comparación del tamaño del genoma mostró que el genoma de O. sinensis era casi 3,4 veces mayor que el de otros hongos entomopatógenos de la familia Hypocreales (Tabla Suplementaria S17 y Figura Suplementaria S8A). El análisis de la secuencia de repetición reveló que esta expansión se debió principalmente a una rápida proliferación de elementos transponibles. Aproximadamente el 74,67% del ensamblaje del genoma de O. sinensis estaba compuesto por secuencias repetidas (Tablas Suplementarias S13-14), excepcionalmente más grandes que las reportadas en Metarhizium anisopliae (~0,98%)19, Metarhizium acridum (~1,52%)19, Cordyceps militaris (~3,04%)20 y Beauveria bassiana (~2,03%)21 (P < 4,822e-07) (Figura Suplementaria S8B). Los elementos MULE fueron notablemente los más abundantes, representando ~1,6% (~1,9 Mb) del genoma de O. sinensis y más del 59% de los transposones de ADN en O. sinensis. Los retrotransposones, en su mayoría retrotransposones de repetición terminal larga (LTR), comprendían el ~59,76% del genoma de O. sinensis y su proliferación a gran escala se produjo aproximadamente a los ~38 millones de años (MYA) (Figura suplementaria S9).

En contraste con la rápida amplificación de los retrotransposones LTR que impulsan la expansión del genoma de O. sinensis, otra característica notable es la dramática pérdida de genes codificadores de proteínas en el linaje de O. sinensis en comparación con otros hongos entomopatógenos. En comparación con un total de 7.939 genes codificadores de proteínas en O. sinensis, había más de 10.095 genes de media en otros hongos entomopatógenos, por ejemplo, Metarhizium anisopliae (10.582)19, Metarhizium acridum (9.849)19, Cordyceps militaris (9.684)20, Beauveria bassiana (10.366)21 y Tolypocladium inflatum (9.998)22 (Tabla 1). Esta reducción del número de genes se evidenció además mediante la identificación de genes no colineales y el análisis comparativo de bloques de sintenia entre los genomas de O. sinensis y C. militaris. Identificamos un total de 308 bloques de sintenia que abarcan casi el 72,7% (~23,4 Mb en C. militaris frente a ~43,5 Mb en O. sinensis) del genoma de C. militaris (Fig. 1A; Figura suplementaria S10 y Tablas suplementarias S18-19). De estas regiones genómicas sintéticas, hubo una disminución de genes no colineales en O. sinensis (2.127) en comparación con C. militaris (3.259) pero un aumento de secuencias repetidas (23,8 Mb en O. sinensis frente a 0,40 Mb en C. militaris) (Fig. 1B y Tabla Suplementaria S19). La anotación funcional de los 2.468 genes que se perdieron en O. sinensis mostró que estaban principalmente implicados en el metabolismo de los aminoácidos, como la biosíntesis de aminoácidos (ko01230), el metabolismo de la arginina y la prolina (ko00330) y el metabolismo de la tirosina (ko00350) (Figura Suplementaria S11 y Tabla Suplementaria S20). En particular, casi el 81% de las secuencias repetidas en estos 308 bloques sintéticos eran retrotransposones LTR, el 40,4% de los cuales eran retroelementos gitanos (Fig. 1B y Tabla Suplementaria S19). La datación molecular estimó que esta clase particular de retrotransposones LTR se amplificó ~38 Mya, lo cual es consistente con el levantamiento de la meseta Qinghai-Tibetana (Fig. 1C).

Evolución rápida de las familias de genes relacionadas con la patogenicidad de los hongos

Una de las características más llamativas del genoma de O. sinensis es la escasez de pares de genes altamente homólogos. De los 7.939 genes predichos que codifican proteínas, ningún par compartía >90% de identidad de aminoácidos en las secuencias de codificación y sólo había un par que compartía >80% de identidad de aminoácidos (Fig. 2A y Tabla Suplementaria S21). Esta característica también se observó en el estrechamente relacionado C. militaris y el hongo ectomicorrícico Tuber melanosporum 23. En comparación con otros hongos entomopatógenos como B. bassiana y C. militaris, las familias multigénicas en O. sinensis eran limitadas en número y comprendían sólo el 8,7% del proteoma predicho; la mayoría de las familias de genes tenían sólo dos miembros (Figura Suplementaria S12). La tasa de ganancia de genes fue sorprendentemente menor que la de pérdida de genes, y entre las 7.800 familias de genes encontradas en el ancestro común más reciente (MRCA) de Hypocreales, 1.756 se perdieron aparentemente en O. sinensis (Fig. 2B). Un espacio de codificación génica tan compacto del genoma de O. sinensis sugiere la naturaleza de este hongo altamente especializado y con poca capacidad de adaptación a múltiples señales ambientales.

Para entender la evolución de las familias de genes que se relacionan con la patogenicidad de los hongos y la adaptación de las tierras altas a entornos difíciles, investigamos las propiedades funcionales de las familias de genes que han sufrido expansiones o contracciones en O. sinensis. El genoma de O. sinensis mostró una considerable expansión de las familias de genes que están principalmente implicadas en la patogenicidad de los hongos, incluyendo la actividad peroxidasa (PF01328; P < 0.01), serina hidrolasa (PF03959; P < 0,01), deuterolisina metaloproteasa (M35) peptidasa (PF02102; P < 0,01), y citocromo P450 (PF00067; P < 0,01) (Tabla Suplementaria S23). Curiosamente, encontramos que las familias de genes ampliadas también están funcionalmente enriquecidas en la categoría Pfam de glucosa-metanol-colina (GMC) oxidorreductasa implicada en el metabolismo ecdysteroid de la muda en los insectos (Tabla Suplementaria S23). En las comparaciones con otros hongos entomopatógenos, la expansión de la familia de genes en el linaje de O. sinensis también se observó con la sobrerrepresentación de términos Pfam putativamente relacionados con la adaptación de la baja temperatura (PF06772; P < 0,01) (Tabla Suplementaria S23).

En cambio, las familias de genes que mostraban un estado de contracción estaban principalmente implicadas en el proceso de transporte y el metabolismo energético, como los transportadores ABC (PF00005; P < 0,01), la permeasa de aminoácidos (PF00324; P < 0,01) y la ATP sintasa (PF00306; P < 0,05) (Tabla Suplementaria S24). Aparte de la evolución dinámica de estas familias de genes, detectamos además 1.077 (~13,57%) genes específicos de la especie en O. sinensis (Fig. 2C). De ellos, 318 (~29,53%) genes pudieron ser anotados funcionalmente y fueron significativamente enriquecidos en categorías GO asociadas con la unión del almidón (GO: 2001070; P < 0,01), la patogénesis (GO: 0009405; P < 0.01), y proceso catabólico de macromoléculas de la pared celular (GO: 0016998; P < 0,01) (Tabla Suplementaria S25).

Para evitar la infección de patógenos fúngicos, los insectos hospedadores suelen producir rápidamente gran cantidad de especies reactivas de oxígeno (ROS) para matar directamente a los patógenos. Como respuesta, los patógenos desarrollaron el sistema de defensa antioxidante de ROS durante la evolución, del cual las peroxidasas, que actúan como enzimas que eliminan ROS, se consideran uno de los componentes más prominentes e integrales24, 25. Entre los genes expandidos en O. sinensis, la actividad peroxidasa fue una de las categorías funcionales altamente enriquecidas (Tabla Suplementaria S23). Las búsquedas en el modelo de Markov oculto (HMM) revelaron 42 (~0,53%) genes de peroxidasa en O. sinensis, cuyo número era notablemente mayor que el de C. militaris (28) y el de la levadura (21), lo que sugiere que una expansión doble de los genes de peroxidasa podría resultar potencialmente en una fuerte capacidad para ayudar a la desintoxicación de ROS en O. sinensis (Fig. 3A y Tabla Suplementaria S26). Entre estos 42 genes de peroxidasas, la haloperoxidasa (haem) es la más abundante, representando ~16,67% del total de peroxidasas detectadas (Fig. 3B). A diferencia de otras especies fúngicas estrechamente relacionadas que carecen por completo de la típica peroxiredoxina de 2-Cisteína, O. sinensis aún conserva una copia (Fig. 3B). La peroxiredoxina 2-Cys se demostró previamente que juega un papel en la respuesta a diferentes niveles de estrés oxidativo en Vibrio vulnificus 26. Un análisis comparativo reveló que el gen retenido en O. sinensis pertenece al tipo Prx1, del que se informó que está funcionalmente conservado27 y que se expresa sólo cuando las células están expuestas a altos niveles de H2O2 generados exógenamente26.

A diferencia del mecanismo de infección de los patógenos de las plantas (PP), que requiere enzimas activas en los carbohidratos (CAZymes) para degradar la compleja pared celular de las plantas28, los patógenos de los insectos (IP) suelen infectar a sus huéspedes penetrando la cutícula29. Para comprobarlo, comparamos O. sinensis y otros cuatro insectos patógenos (M. anisopliae, M. acridum, C. militaris y B. bassiana) con los cuatro patógenos de plantas (Fusarium graminearum, Magnaporthe grisea, Grosmannia clavigera y Botrytis cinerea) (Tabla suplementaria S17). Nuestros resultados demostraron que los patógenos de insectos tenían más proteasas (una media de 362 en los PI frente a 342 en los PP; P < 0,43) y proteínas quinasas (una media de 151 en los PI frente a 119 en los PP; P < 0,0014) para degradar la cutícula de los insectos en comparación con los patógenos de las plantas (Tablas suplementarias S27-29). Por el contrario, los patógenos de las plantas albergaban más CAZimas que los patógenos de los insectos para la degradación de la pared celular de las plantas (en promedio, 161 en los PI frente a 231 en los PP) (Tablas Suplementarias S30-32). Excluyendo los patógenos de plantas, O. sinensis tenía notablemente menos genes que codifican proteasas (260) que otros patógenos de insectos, como M. anisopliae (437), M. acridum (361), C. militaris (355) y B. bassiana (396). Sin embargo, el ~35% de estas proteasas caracterizadas en O. sinensis contienen un péptido señal que es más probable que haya participado en las interacciones patógeno-hospedador (Tablas Suplementarias S10 y S34), lo que es mayor que en otros hongos entomopatógenos (una media del 20%). Al igual que en otros patógenos de insectos, varias familias de celulasa, incluyendo GH7, GH45 y GH51, también disminuyeron o estuvieron ausentes en O. sinensis (Tabla Suplementaria S30).

También examinamos los perfiles de expresión génica en los tres estadios de desarrollo de O. sinensis con ratios de longitud del hongo frente al insecto que alcanzaron ~1,20×, ~1,75× y ~2,20×. Los resultados muestran que un total de 411 genes fueron expresados diferencialmente (DEG) entre las tres etapas de desarrollo (Figura Suplementaria S14). La anotación funcional de estos 411 DEGs descubrió que estaban principalmente implicados en la patogenicidad de los hongos, como la familia 16 de glicosil hidrolasas (PF00722; FDR < 0,01), el citocromo P450 (PF00067; FDR < 0,01) y la superfamilia de facilitadores principales (PF07690; FDR < 0,05). Además, los genes que codifican enzimas asociadas a la cadena respiratoria mitocondrial también estaban funcionalmente enriquecidos, como la familia de la epimerasa/deshidratasa dependiente de NAD (PF01370; FDR < 0,01) y el dominio N-terminal de BCS1 (PF08740; FDR < 0,01) (Tabla suplementaria S33).

La selección darwiniana positiva sirve como fuerza motriz para la patogenicidad de los hongos

La selección positiva ha desempeñado sin duda un papel crítico en la evolución de diversos organismos que viven en entornos de gran altitud de la meseta Qinghai-Tibetano, y es probable que muchos de los rasgos fenotípicos muestren tales firmas de selección3,4,5. De los 1.499 ortólogos de alta confianza compartidos entre O. sinensis y las otras 12 especies de hongos, se identificaron 163 genes positivamente seleccionados (PSGs) en O. sinensis utilizando la prueba de relación de verosimilitud de la rama (LRT; P < 0,05) (Tabla Suplementaria S35). De ellos, un gen (OSIN3929; aquí denominado OsPRX1) ha sido implicado funcionalmente en la actividad de la peroxidasa (Fig. 3C). Este gen es un miembro de la familia de las peroxiredoxinas con 1-cisteína y altamente homólogo a PRX1 (YBL064C) en S. cerevisiae 30. Los genes PRX1 en S. cerevisiae y en dos hongos patógenos humanos, A. fumigatus y C. albicans, están funcionalmente conservados y son necesarios para la desintoxicación de la explosión oxidativa dentro de las células huésped31, 32. En particular, la supresión de PRX1 en el conocido patógeno del arroz, Magnaporthe oryzae, dio lugar a una pérdida casi completa de patogenicidad, lo que sugiere que esta peroxidasa es clave para las interacciones huésped-patógeno27. Sorprendentemente, también se detectó que varios genes implicados en las interacciones huésped-patógeno, incluyendo la biogénesis peroxisomal, la proteína quinasa y las metalopéptidasas, estaban bajo selección positiva (Fig. 3C).

Evolución del sistema de apareamiento

En los hongos ascomicetos, el sistema de apareamiento suele estar controlado por el locus de apareamiento (MAT)33. Nuestro análisis de secuenciación del genoma encontró que O. sinensis no sólo poseía el gen de tipo de apareamiento MAT1-2-1 dentro del idiomorfo MAT1-2, sino que también tenía tres genes de tipo de apareamiento (es decir, MAT1-1-1, MAT1-1-2 y MAT1-1-3) dentro del idiomorfo MAT1-1 (Figura Suplementaria S15B). Esta característica se verificó utilizando la resecuenciación del genoma completo de 31 poblaciones naturales en casi toda el área de distribución geográfica, lo que indica que O. sinensis es realmente homotálica (Figura S15A suplementaria y Tabla S36 suplementaria). Esta característica es extremadamente diferente de sus hongos patógenos estrechamente relacionados, como Tolypocladium inflatum (MAT1-2)22, C. militaris (MAT1-1)20, B. bassiana (MAT1-1)21, M. anisopliae (MAT1-1)19, y M. acridum (MAT1-2)19, que son heterotálicos y poseen un solo locus de apareamiento. Al igual que el conocido patógeno vegetal homotálico Fusarium graminearum 34, la organización de estos dos loci MAT en O. sinensis reveló el estado de fusión dentro de la región genómica idiomorfa, que estaba particularmente enriquecida en retrotransponsones LTR. La escisión entre el O. sinensis homotálico y el C. militaris heterotálico se estimó que ocurrió cerca de 174,2 MYA (Figura Suplementaria S13C), y fue objeto de múltiples conversiones del sistema de apareamiento de autoincompatible a autocompatible durante el curso de su historia evolutiva (Figura Suplementaria S15C), asemejándose al género de ascomicetos filamentosos de Neurospora 35.

Divergencia poblacional basada en las latitudes de la meseta Qinghai-Tibetana

Para examinar las relaciones en todo el genoma y la divergencia poblacional, recogimos y resecuenciamos 31 accesiones de O. sinensis a lo largo de su rango de distribución conocido, incluyendo las provincias de Qinghai, Sichuan, Yunnan y Gansu y la Región Autónoma del Tíbet en la meseta Qinghai-Tibetana (Figura Suplementaria S16 y Tabla Suplementaria S37). Generamos un total de 183 millones de lecturas pareadas (~36,68 Gb de secuencias) con una profundidad media de ~10,1× (datos brutos) (Tabla Suplementaria S38). A partir de estos datos, generamos un conjunto de 816.960 polimorfismos de un solo nucleótido (SNP) y 48.092 indels estrictos (inserciones y deleciones) para evaluar el parentesco entre las poblaciones de O. sinensis (Figuras suplementarias S18-19 y Tabla suplementaria S39). La mayoría de las variantes genómicas (71,1%) se ubicaron en las regiones intergénicas, con un subconjunto que se ubicó en las regiones codificantes (23,3%, consistente en 101.997 SNP sinónimos y 88.224 no sinónimos, con una proporción de sustitución de 0,86) (Figura Suplementaria S19 y Tabla Suplementaria S39). El árbol filogenético construido utilizando los conjuntos de datos SNP dividió las 31 accesiones en tres grupos geográficamente separados que van desde las regiones de baja latitud a las de alta latitud (Fig. 4A), un hallazgo que fue reforzado por el PCA utilizando los primeros y segundos vectores propios (Fig. 4B y Figura Suplementaria S20A). La variación del número de presuntas poblaciones ancestrales (K) mostró que cuando K = 3, los tres grupos distintos son consistentes con el PCA y la reconstrucción filogenética (Fig. 4C y Figura Suplementaria S20B). Algunas accesiones del grupo de baja latitud muestran una fuerte evidencia de la mezcla y están más dispersas en comparación con los otros dos grupos, indicando una mayor diversidad genética, posiblemente debido a los polimorfismos ancestrales compartidos y/o eventos de introgresión recientes (Fig. 4D,E). La estadística de diferenciación poblacional estimada (F ST ) entre estos tres grupos basados en la latitud reveló aún más la naturaleza basal de la región de baja latitud, particularmente de las poblaciones del distrito de Nyingchi en el Tíbet, lo que se evidenció además por su diversidad de nucleótidos sustancialmente elevada (π) dentro del grupo y una menor diferenciación poblacional con los otros dos grupos de alta latitud. (Fig. 4C-F).

Investigamos además los genes afectados por diferentes niveles de contenidos de SNP y mutaciones no sinónimas (Tabla Suplementaria S40). El análisis de enriquecimiento funcional de los 100 genes principales con el mayor contenido de SNP y/o mutaciones no sinónimas muestra que están principalmente implicados en el metabolismo de los metabolitos secundarios de los hongos, como la sintasa deshidratada de policétidos (PF14765; FDR < 0,01), el dominio KR (PF08659; FDR < 0,01) y el dominio de condensación (PF00668; FDR < 0,01). También se enriquecieron las categorías funcionales asociadas a la biosíntesis de ácidos grasos, como el dominio de la acil transferasa (PF00698; FDR < 0,01) y la beta-cetoacil sintasa (PF00109 y PF02801; FDR < 0,01) (Tablas suplementarias S41-42).