Séquençage du génome, assemblage et annotation

Nous avons séquencé O. sinensis du district de Nyingchi au Tibet, en Chine. Nous avons effectué une analyse WGS avec les plateformes de séquençage de nouvelle génération Roche 454 et Illumina HiSeq 2000. Cela a généré des ensembles de données de séquences propres de ~5,4 Gb, ce qui donne une couverture du génome d’environ 45,1 fois, respectivement (Tableau supplémentaire S1). Nous avons estimé que la taille du génome est de ~124,08 Mb et ~119,8 Mb sur la base de la cytométrie en flux et de la distribution de la profondeur 17-mer des lectures séquencées, respectivement (Figures supplémentaires S1-2 et Tableau supplémentaire S2). Le génome de O. sinensis a d’abord été assemblé à partir des longues lectures de Roche 454 à l’aide de Newbler16, puis les contigs pré-assemblés ont été échafaudés avec les lectures de séquençage de paires appariées d’Illumina à l’aide de SSPACE17. Cela a finalement donné un assemblage du génome de ~116,4 Mb qui couvre ~97 % de la taille estimée du génome et contient 156 échafaudages (>2 Kb) avec une valeur ScafN50 de ~3 Mb et 9 141 contigs (N50 = 21 423 pb) (tableau 1 et tableaux supplémentaires S3-4). Pour valider la qualité de l’assemblage du génome, nous avons d’abord aligné tout l’ADN et les étiquettes de séquences exprimées (EST) de O. sinensis disponibles dans les bases de données publiques et avons obtenu des taux de cartographie de 98,85 % et 95,33 %, respectivement (Tableau supplémentaire S5). Deuxièmement, nous avons mappé toutes les longues lectures Roche 454 propres (~1,84 Gb) aux séquences du génome assemblé, montrant un alignement presque parfait avec un taux de mappage de 99,01% (Tableau supplémentaire S5). Troisièmement, les transcriptions que nous avons assemblées ont montré un bon alignement sur le génome assemblé ; sur 11 742 transcriptions, 91,29 % ont été cartographiées (couverture des transcriptions ≥80 % et identité ≥90 % ; tableau supplémentaire S5). Enfin, nous avons évalué la complétude de notre assemblage O. sinensis en utilisant BUSCO18 ; 94,0 %, 4,0 % et 1,8 % des 1 315 gènes conservés BUSCO Ascomycota attendus ont été identifiés comme complets, fragmentés et manquants dans notre assemblage O. sinensis, respectivement (tableau supplémentaire S5).

Nous avons généré ~15,05 Go de données de séquençage de l’ARN (RNA-Seq) obtenues à partir d’un total de six bibliothèques représentant les trois principaux stades de développement pour faciliter la prédiction des gènes (figure supplémentaire S4 et tableaux supplémentaires S6,7). En combinaison avec la prédiction ab initio, les alignements de protéines et d’EST, le peignage EvidenceModeler et un filtrage supplémentaire, nous avons défini 7 939 gènes codant pour des protéines (tableau 1 et tableau supplémentaire S8). Parmi ces gènes prédits, environ 97,0 % et 71,51 % ont pu être classés fonctionnellement et étayés par des données RNA-Seq, respectivement (tableaux supplémentaires S9-11). En utilisant le BUSCO de la lignée Ascomycota, nous avons également constaté que 94,4 %, 3,6 %, 1,8 % et 0,2 % des gènes étaient complets, fragmentés, manquants et dupliqués, respectivement, ce qui indique une bonne qualité de notre annotation génétique (tableau supplémentaire S11). Nous avons également effectué des recherches d’homologues et annoté les gènes d’ARN non codant (ARNnc), ce qui nous a permis d’identifier 146 gènes d’ARN de transfert (ARNt), 33 gènes d’ARN ribosomal (ARNr), 70 gènes de petits ARN nucléolaires (snoRNA) et 15 gènes de petits ARN nucléaires (snRNA) (figure supplémentaire S6 et tableau supplémentaire S12). L’annotation des séquences répétées a montré que les éléments transposables (TE) représentaient environ 74,67 % du génome assemblé et 80,07 % des lectures brutes, ce qui indique que ~5,45 % du génome non assemblé est constitué de TE (tableaux supplémentaires S13-14). La teneur en GC était de 43,09 % dans l’ensemble du génome et de 61,49 % dans les séquences codantes (figure supplémentaire S3 ; tableaux supplémentaires S4 et S8). Nous avons annoté 8 918 répétitions de séquences simples qui fourniront des marqueurs génétiques précieux pour aider les futurs programmes de sélection du champignon de la chenille chinoise (tableaux supplémentaires S15-16 et figure supplémentaire S7).

Expansion du génome induite par les rétrotransposons et élimination massive des gènes non colinéaires

La comparaison de la taille des génomes a montré que le génome d’O. sinensis était près de 3,4 fois plus grand que celui des autres champignons entomopathogènes de la famille des Hypocreales (tableau supplémentaire S17 et figure supplémentaire S8A). L’analyse des séquences répétées a révélé que cette expansion était principalement due à une prolifération rapide d’éléments transposables. Environ 74,67 % de l’assemblage du génome d’O. sinensis était composé de séquences répétées (tableaux supplémentaires S13-14), ce qui est exceptionnellement plus important que ceux signalés chez Metarhizium anisopliae (~0,98 %)19, Metarhizium acridum (~1,52 %)19, Cordyceps militaris (~3,04 %)20 et Beauveria bassiana (~2,03 %)21 (P < 4,822e-07) (figure supplémentaire S8B). Les éléments MULE étaient notamment les plus abondants, représentant ~1,6% (~1,9 Mb) du génome d’O. sinensis et plus de 59% des transposons d’ADN chez O. sinensis. Les rétrotransposons, principalement les rétrotransposons à longue répétition terminale (LTR), représentaient ~59,76% du génome de O. sinensis et dont la prolifération à grande échelle s’est produite approximativement à ~38 millions d’années (MYA) (Supplemental Figure S9).

Évolution rapide des familles de gènes liées à la pathogénicité fongique

L’une des caractéristiques les plus frappantes du génome d’O. sinensis est la pénurie de paires de gènes hautement homologues. Parmi les 7 939 gènes codant pour des protéines, aucune paire ne partageait >90% d’identité d’acides aminés dans les séquences codantes et une seule paire partageait >80% d’identité d’acides aminés (Fig. 2A et Tableau supplémentaire S21). Cette caractéristique a également été observée chez le C. militaris étroitement apparenté et le champignon ectomycorhizien Tuber melanosporum 23. Par rapport à d’autres champignons entomopathogènes comme B. bassiana et C. militaris, les familles multigéniques chez O. sinensis étaient limitées en nombre et ne représentaient que 8,7 % du protéome prédit ; la plupart des familles de gènes ne comptaient que deux membres (figure supplémentaire S12). Le taux de gain de gènes était étonnamment plus faible que celui de perte de gènes, et parmi les 7 800 familles de gènes trouvées dans l’ancêtre commun le plus récent (MRCA) des Hypocreales, 1 756 étaient apparemment perdues chez O. sinensis (Fig. 2B). Un tel espace compact de codage des gènes du génome d’O. sinensis suggère la nature de ce champignon hautement spécialisé avec une faible capacité d’adaptation à de multiples signaux environnementaux.

Pour comprendre l’évolution des familles de gènes qui sont liées à la pathogénicité fongique et à l’adaptation des hautes terres à des environnements difficiles, nous avons étudié les propriétés fonctionnelles des familles de gènes qui ont subi des expansions ou des contractions chez O. sinensis. Le génome de O. sinensis a montré une expansion considérable des familles de gènes qui sont principalement impliqués dans la pathogénicité fongique, y compris l’activité peroxydase (PF01328 ; P < 0.01), la sérine hydrolase (PF03959 ; P < 0,01), la deutérolysine métalloprotéase (M35) peptidase (PF02102 ; P < 0,01), et le cytochrome P450 (PF00067 ; P < 0,01) (Tableau supplémentaire S23). De manière intéressante, nous avons constaté que les familles de gènes élargies sont également enrichies fonctionnellement dans la catégorie Pfam de la glucose-méthanol-choline (GMC) oxydoréductase impliquée dans le métabolisme ecdystéroïde de la mue chez les insectes (Tableau supplémentaire S23). En comparant avec d’autres champignons entomopathogènes, l’expansion de la famille de gènes dans la lignée de O. sinensis a également été observée avec une surreprésentation des termes Pfam putativement liés à l’adaptation de la basse température (PF06772 ; P < 0,01) (Tableau supplémentaire S23).

Pour éviter l’infection des champignons pathogènes, les insectes hôtes produisent souvent rapidement beaucoup d’espèces réactives de l’oxygène (ROS) pour tuer directement les pathogènes. En réponse, les agents pathogènes ont évolué le système de défense antioxydant ROS au cours de l’évolution, dont les peroxydases, agissant comme des enzymes de piégeage ROS, sont considérés comme l’un des composants les plus importants et intégral24, 25. Parmi les gènes développés chez O. sinensis, l’activité peroxydase était l’une des catégories fonctionnelles hautement enrichies (tableau supplémentaire S23). Les recherches par modèle de Markov caché (HMM) ont révélé 42 (~0,53%) gènes de peroxydase chez O. sinensis, dont le nombre était remarquablement plus élevé que celui de C. militaris (28) et de la levure (21), ce qui suggère qu’une multiplication par deux des gènes de peroxydase pourrait potentiellement se traduire par une forte capacité à contribuer à la détoxification des ROS chez O. sinensis (Fig. 3A et Tableau supplémentaire S26). Parmi ces 42 gènes de peroxydase, l’haloperoxydase (hème) est la plus abondante, représentant environ 16,67 % du total des peroxydases détectées (Fig. 3B). Contrairement à d’autres espèces fongiques étroitement apparentées qui manquent complètement de la peroxiredoxine 2-Cystéine typique, O. sinensis conserve encore une copie (Fig. 3B). Il a été démontré précédemment que la peroxiredoxine 2-Cys joue un rôle dans la réponse à différents niveaux de stress oxydatif chez Vibrio vulnificus 26. Une analyse comparative a révélé que le gène retenu chez O. sinensis appartient au type Prx1, qui a été signalé comme étant fonctionnellement conservé27 et exprimé uniquement lorsque les cellules sont exposées à des niveaux élevés de H2O2 générés de manière exogène26.

Contrairement au mécanisme d’infection des agents pathogènes des plantes (PPs), qui nécessite des enzymes à hydrate de carbone-actif (CAZymes) pour dégrader la paroi cellulaire complexe des plantes28, les insectes pathogènes (IPs) infectent généralement leurs hôtes en pénétrant la cuticule29. Pour tester cela, nous avons comparé O. sinensis et quatre autres insectes pathogènes (M. anisopliae, M. acridum, C. militaris et B. bassiana) avec les quatre agents pathogènes des plantes (Fusarium graminearum, Magnaporthe grisea, Grosmannia clavigera et Botrytis cinerea) (Tableau supplémentaire S17). Nos résultats ont démontré que les agents pathogènes des insectes possédaient davantage de protéases (en moyenne, 362 dans les IP contre 342 dans les PP ; P < 0,43) et de protéines kinases (en moyenne, 151 dans les IP contre 119 dans les PP ; P < 0,0014) pour dégrader la cuticule des insectes par rapport aux agents pathogènes des plantes (tableaux supplémentaires S27-29). En revanche, les agents pathogènes des plantes abritaient plus de CAZymes que les insectes pathogènes pour la dégradation de la paroi cellulaire des plantes (en moyenne, 161 dans les IP contre 231 dans les PP) (tableaux supplémentaires S30-32). Si l’on exclut les agents pathogènes végétaux, O. sinensis possède remarquablement moins de gènes codant pour des protéases (260) que d’autres insectes pathogènes, tels que M. anisopliae (437), M. acridum (361), C. militaris (355) et B. bassiana (396). Cependant, ~35% de ces protéases caractérisées chez O. sinensis contiennent un peptide signal qui est plus susceptible d’avoir été impliqué dans les interactions pathogène-hôte (Tableaux supplémentaires S10 et S34), ce qui est plus important que chez les autres champignons entomopathogènes (en moyenne, 20%). Comme pour les autres insectes pathogènes, plusieurs familles de cellulases, dont GH7, GH45 et GH51, ont également diminué ou étaient absentes chez O. sinensis (tableau supplémentaire S30).

Nous avons également examiné les profils d’expression génique à travers les trois stades de développement d’O. sinensis, les rapports de longueur du champignon par rapport à l’insecte atteignant ~1,20×, ~1,75× et ~2,20×. Les résultats montrent qu’un total de 411 gènes ont été exprimés de manière différentielle (DEG) entre les trois stades de développement (figure supplémentaire S14). L’annotation fonctionnelle de ces 411 DEG a révélé qu’ils étaient principalement impliqués dans la pathogénicité fongique, comme les glycosyl hydrolases de la famille 16 (PF00722 ; FDR < 0,01), le cytochrome P450 (PF00067 ; FDR < 0,01) et la superfamille des facilitateurs majeurs (PF07690 ; FDR < 0,05). En outre, les gènes codant pour les enzymes associées à la chaîne respiratoire mitochondriale étaient également enrichis sur le plan fonctionnel, comme la famille des épimérases/déshydratases dépendantes du NAD (PF01370 ; FDR < 0,01) et le domaine N-terminal de BCS1 (PF08740 ; FDR < 0,01) (Tableau supplémentaire S33).

La sélection darwinienne positive sert de forces motrices pour la pathogénicité fongique

La sélection positive a sans aucun doute joué un rôle critique dans l’évolution de divers organismes vivant dans les environnements de haute altitude du plateau Qinghai-Tibetan, et de nombreux traits phénotypiques sont susceptibles de présenter de telles signatures de sélection3,4,5. Sur les 1 499 orthologues monocopies de haute confiance partagés entre O. sinensis et les 12 autres espèces fongiques, 163 gènes sélectionnés positivement (PSG) ont été identifiés chez O. sinensis en utilisant le test du rapport de vraisemblance (LRT ; P < 0,05) (Tableau supplémentaire S35). Parmi eux, un gène (OSIN3929 ; ici nommé OsPRX1) a été impliqué fonctionnellement dans l’activité peroxydase (Fig. 3C). Ce gène est un membre de la famille des peroxiredoxines avec 1-cystéine et hautement homologue à PRX1 (YBL064C) dans S. cerevisiae 30. Les gènes PRX1 de S. cerevisiae et de deux champignons pathogènes pour l’homme, A. fumigatus et C. albicans, sont conservés sur le plan fonctionnel et nécessaires à la détoxification de la salve oxydative dans les cellules hôtes31, 32. En particulier, la délétion de PRX1 dans l’agent pathogène bien connu du riz, Magnaporthe oryzae, a entraîné une perte presque complète de pathogénicité, ce qui suggère que cette peroxydase est essentielle aux interactions hôte-pathogène27. De façon frappante, plusieurs gènes impliqués dans les interactions hôte-pathogène, y compris la biogenèse peroxysomale, la protéine kinase et les métallopeptidases, ont également été détectés pour être sous sélection positive (Fig. 3C).

Évolution du système d’accouplement

Dans les champignons ascomycètes, le système d’accouplement est généralement contrôlé par le locus de type d’accouplement (MAT)33. Notre analyse de séquençage du génome a révélé que O. sinensis possédait non seulement le gène de type d’accouplement MAT1-2-1 au sein de l’idiomorphe MAT1-2, mais aussi trois gènes de type d’accouplement (c’est-à-dire MAT1-1-1, MAT1-1-2 et MAT1-1-3) au sein de l’idiomorphe MAT1-1 (figure supplémentaire S15B). Cette caractéristique a été vérifiée à l’aide du reséquençage du génome entier de 31 populations naturelles dans la quasi-totalité de l’aire géographique, ce qui indique que O. sinensis est effectivement homothallique (Supplemental Figure S15A et Supplemental Table S36). Cette caractéristique est extrêmement différente de celle des champignons pathogènes étroitement apparentés, tels que Tolypocladium inflatum (MAT1-2)22, C. militaris (MAT1-1)20, B. bassiana (MAT1-1)21, M. anisopliae (MAT1-1)19 et M. acridum (MAT1-2)19, qui sont hétérogalliques et ne possèdent qu’un seul locus de type d’accouplement. Comme pour l’agent phytopathogène homothallique bien connu Fusarium graminearum 34, l’organisation de ces deux loci MAT chez O. sinensis a révélé l’état de fusion au sein de la région génomique idiomorphe, qui était particulièrement enrichie en rétrotransponsons LTR. La scission entre l’homothallique O. sinensis et l’hétérothallique C. militaris a été estimée à près de 174,2 MYA (figure supplémentaire S13C), et a été soumise à de multiples conversions du système d’accouplement d’auto-incompatible à auto-compatible au cours de son histoire évolutive (figure supplémentaire S15C), ressemblant au genre ascomycète filamenteux de Neurospora 35.

Divergence de la population en fonction des latitudes sur le plateau Qinghai-Tibétain

Pour examiner les relations à l’échelle du génome et la divergence de la population, nous avons recueilli et reséquencé 31 accessions d’O. sinensis dans toute son aire de répartition connue, y compris les provinces de Qinghai, Sichuan, Yunnan et Gansu et la région autonome du Tibet sur le plateau Qinghai-Tibétain (figure supplémentaire S16 et tableau supplémentaire S37). Nous avons généré un total de 183 millions de lectures en paires (~36,68 Go de séquences) avec une profondeur moyenne de ~10,1× (données brutes) (Tableau supplémentaire S38). À partir de ces données, nous avons généré un ensemble de 816 960 polymorphismes mononucléotidiques (SNP) et 48 092 indels stricts (insertions et délétions) pour évaluer la parenté entre les populations d’O. sinensis (Figures supplémentaires S18-19 et Tableau supplémentaire S39). La majorité des variants génomiques (71,1 %) ont été cartographiés dans les régions intergéniques avec un sous-ensemble cartographié dans les régions codantes (23,3 % consistant en 101 997 SNP synonymes et 88 224 SNP non synonymes avec un taux de substitution de 0,86) (figure supplémentaire S19 et tableau supplémentaire S39). L’arbre phylogénétique construit à l’aide des ensembles de données SNP a divisé les 31 entrées en trois groupes géographiquement distincts allant des régions de basse latitude aux régions de haute latitude (Fig. 4A) – une conclusion qui a été renforcée par l’ACP utilisant les premier et deuxième vecteurs propres (Fig. 4B et Supplemental Figure S20A). La variation du nombre de populations ancestrales présumées (K) a montré que lorsque K = 3, les trois groupes distincts sont cohérents avec l’ACP et la reconstruction phylogénétique (Fig. 4C et Supplemental Figure S20B). Certaines entrées du groupe des basses latitudes présentent de fortes preuves de mélange et sont plus dispersées que celles des deux autres groupes, ce qui indique une plus grande diversité génétique probablement due à des polymorphismes ancestraux partagés et/ou à des événements d’introgression récents (Fig. 4D,E). La statistique de différenciation de la population estimée (F ST ) parmi ces trois groupes basés sur la latitude a en outre révélé la nature basale de la région des basses latitudes, en particulier les populations du district de Nyingchi au Tibet, ce qui a été mis en évidence par sa diversité nucléotidique sensiblement élevée (π) au sein du groupe et une différenciation de la population plus faible avec les deux autres groupes des hautes latitudes. (Fig. 4C-F).

Nous avons ensuite étudié les gènes affectés par différents niveaux de contenu SNP et de mutations non synonymes (tableau supplémentaire S40). L’analyse d’enrichissement fonctionnel des 100 principaux gènes présentant le contenu SNP et/ou les mutations non synonymes les plus élevés montre qu’ils sont principalement impliqués dans le métabolisme des métabolites secondaires fongiques, tels que la polykétide synthase déshydratase (PF14765 ; FDR < 0,01), le domaine KR (PF08659 ; FDR < 0,01) et le domaine de condensation (PF00668 ; FDR < 0,01). Les catégories fonctionnelles associées à la biosynthèse des acides gras étaient également enrichies, telles que le domaine acyltransférase (PF00698 ; FDR < 0,01) et la bêta-cétoacyl-synthase (PF00109 et PF02801 ; FDR < 0,01) (Tableaux supplémentaires S41-42).