Genome sequencing, assembly and annotation

Nyingchi District of Tibet, ChinaのO. sinensisをシークエンスしました。 次世代シーケンサーRoche 454とIllumina HiSeq 2000のプラットフォームでWGS解析を行った。 この結果、約5.4Gbのクリーンな配列データセットが得られ、それぞれ約45.1倍のゲノムカバレッジが得られた（Supplemental Table S1）。 フローサイトメトリーおよび配列決定されたリードの17-mer深さ分布から、ゲノムサイズはそれぞれ〜124.08 Mbおよび〜119.8 Mbと推定されました（補足図S1-2および補足表S2）。 O. sinensisゲノムは、まずNewbler16を用いてRoche 454ロングリードからアセンブルし、次にSSPACE17を用いてイルミナメイトペアシーケンスリードで事前にアセンブルしたコンティグを足場としました。 この結果、最終的に推定ゲノムサイズの97%をカバーする約116.4Mbのゲノムアセンブリが得られ、ScafN50値が約3Mbの156個のスキャフォールド（>2Kb）と9141個のコンティグ（N50 = 21,423 bp）が含まれていた（表1および補足表S3-4）。 ゲノムアセンブリの品質を検証するために、まず、公共データベースで利用可能なO. sinensisのすべてのDNAおよび発現配列タグ（EST）をアラインメントし、それぞれ98.85%および95.33%のマッピング率を得た（Supplemental Table S5）。 次に、Roche 454ロングリード（約1.84 Gb）をすべてアセンブルしたゲノム配列にマッピングしたところ、ほぼ完璧なアライメントを示し、マッピング率は99.01%でした（Supplemental Table S5）。 第三に、我々が組み立てた転写産物は、組み立てられたゲノムに対して良好なアライメントを示し、11,742の転写産物のうち91.29%がマッピングされました（転写産物カバレッジ≧80%、同一性≧90%、Supplemental Table S5)。 最後に、BUSCO18を用いてO. sinensisアセンブリの完全性を評価したところ、1,315の予想される子嚢菌BUSCO保存遺伝子の94.0%、4.0%、1.8%がそれぞれO. sinensisアセンブリにおいて完全、断片、欠損として識別されました（Supplemental Table S5）。

遺伝子予測を助けるために、3つの主要な発生段階を表す合計6ライブラリーから得られた約15.05 GbのRNAシーケンス（RNA-SEQ）データを作成しました（補足図S4および補足表S6,7）。 ab initio予測、タンパク質とESTのアラインメント、EvidenceModelerの結合、さらにフィルタリングを組み合わせて、7,939個のタンパク質コード遺伝子を定義しました（表1、補足表S8）。 これらの予測された遺伝子のうち、約97.0%と71.51%が機能的に分類され、それぞれRNA-Seqデータでサポートされました（補足表S9-11）。 また、子嚢菌類由来のBUSCOを用いた場合、94.4%、3.6%、1.8%、0.2%が完全、断片化、欠損、重複しており、我々の遺伝子アノテーションは良好な品質を持つことが分かりました（Supplemental Table S11）。 また、相同性検索を行い、非コードRNA（ncRNA）遺伝子のアノテーションを行った結果、トランスファーRNA（tRNA）遺伝子146個、リボソームRNA（rRNA）遺伝子33個、小核RNA（snoRNA）遺伝子70個、小核RNA（snRNA）遺伝子15個が得られた（補足図S6および補足表S12）。 繰り返し配列のアノテーションでは、トランスポーザブル・エレメント（TE）が、組み立てられたゲノムの約74.67％、生のリードの80.07％を占め、未組み立てゲノムの約5.45％がTEからなることが示された（補足表S13〜14）。 GC含量はゲノム全体で43.09%、コーディング配列で61.49%でした（補足図S3; 補足表S4, S8）。 また、8,918個の単純配列反復配列をアノテーションし、将来の中国イモムシ菌の育種プログラムを支援するための貴重な遺伝マーカーを提供した（補足表S15〜16、補足図S7）。

レトロトランスポゾンによるゲノム拡大と非共線性遺伝子の大量除去

ゲノムサイズの比較から、オサムシゲノムがヒポクラゲ科の他の昆虫病原性菌に比べて約3.4倍大きいことがわかった（補足表S17、補足図S8a）。 繰り返し配列解析の結果、この拡大は主にトランスポーザブルエレメントの急速な増殖によるものであることが明らかになった。 O. sinensisゲノムアセンブリの約74.67%がリピート配列で構成されており（補足表S13-14）、Metarhizium anisopliae（~0.98%）19, Metarhizium acridum（~1.52%）19, Cordyceps militaris（~3.04%）20 and Beauveria bassiana（~2.03%）21 (P < 4.822e-07) より格段大きくなった(補足図S8B)。 MULEエレメントが最も多く、O. sinensisゲノムの約1.6% (~1.9 Mb)を占め、O. sinensisのDNAトランスポゾンの59%以上を占めたことは特筆に値する。 レトロトランスポゾンは、主にLTRレトロトランスポゾンで、O. sinensisゲノムの〜59.76%を占め、その大規模な増殖はおよそ3800万年前（MYA）に起こった（補足図S9）。

LTRレトロトランスポゾンの急速な増幅がO. sinensisゲノムの拡大を促進したのとは対照的に、O. sinensis系統では他の昆虫病原真菌と比較してタンパク質コード遺伝子が劇的に失われていることも顕著な特徴であった。 例えば、Metarhizium anisopliae (10,582)19, Metarhizium acridum (9,849)19, Cordyceps militaris (9,684)20, Beauveria bassiana (10,366)21 and Tolypocladium inflatum (9,998)22 (Table 1)などの昆虫病原性細菌では平均で 10,095 遺伝子以上存在していた（O. sinensis では合計 7,939 遺伝子のうち、10,095 遺伝子が失われていた）。 このような遺伝子数の減少は、非共線性遺伝子の同定やO. sinensisとC. militarisゲノムのシンテニーブロックの比較解析によってさらに証明された。 その結果、C. militarisゲノムの約72.7%（C. militarisでは約23.4Mb、O. sinensisでは約43.5Mb）に及ぶ308のシンテニックブロックを特定した（図1A；補足図S10および補足表S18-19）。 これらのシンテニックゲノム領域のうち、O. sinensisではC. militaris（3,259）よりも非共線性遺伝子が減少していたが（O. sinensis 23.8 Mb vs. C. militaris 0.40 Mb）、繰り返し配列が増加していた（図1Bおよび補足表S19）。 O. sinensisで失われた2,468個の遺伝子の機能アノテーションによると、それらは主にアミノ酸の生合成（ko01230）、アルギニンおよびプロリン代謝（ko00330）、チロシン代謝（ko00350）などアミノ酸代謝に関与していた（補足図S11、補足表S20）。 注目すべきは、これら308のシンテニックブロックのリピート配列の81％近くがLTRレトロトランスポゾンで、そのうち40.4％がジプシーレトロエレメントだったことである（図1B、補足表S19）。 分子年代測定では、この特定のクラスのLTRレトロトランスポゾンは38Maに増幅されたと推定され、これは青海チベット高原の隆起と一致する（図1C）。

真菌の病原性に関わる遺伝子ファミリーの急速な進化

O. sinensisゲノムで最も顕著な特徴の一つは相同性の高い遺伝子対がほとんど存在しないことである。 予測された7,939個のタンパク質コード遺伝子のうち、コード配列において>90%のアミノ酸同一性を持つペアはなく、>80%のアミノ酸同一性を持つペアは1つしかなかった（図2Aおよび補足表S21）。 この特徴は、近縁種のC. militarisや外生菌根菌のTuber melanosporum 23でも観察された。 B. bassianaやC. militarisのような他の昆虫病原真菌と比較して、O. sinensisの多遺伝子ファミリーは数が少なく、予測されるプロテオームの8.7%を占めるのみで、ほとんどの遺伝子ファミリーは2メンバーのみだった（補遺図S12）。 遺伝子の増加率は減少率より著しく低く、ハイポクレアーレの最新共通祖先（MRCA）で見つかった7,800の遺伝子ファミリーのうち、1,756がO. sinensisでは失われたようである（図2B)。 O. sinensisゲノムのこのようなコンパクトな遺伝子コード空間は、複数の環境キューに適応する能力の低い、この高度に特殊化した菌の性質を示唆している。

真菌の病原性や過酷な環境への高地適応に関わる遺伝子ファミリーの進化を理解するために、O. sinensisで拡大または縮小を遂げた遺伝子ファミリーの機能特性を調査した。 O. sinensisゲノムでは、主に真菌の病原性に関わる遺伝子ファミリーがかなり拡大し、ペルオキシダーゼ活性（PF01328；P < 0.01）、セリンヒドロラーゼ（PF03959；P < 0.01）、デューテロライシンメタロプロテアーゼ（M35）ペプチダーゼ（PF02102；P < 0.01）、シトクロムP450（PF00067；P < 0.01）（Suplemental Table S23）であった。 興味深いことに、拡張された遺伝子ファミリーは、昆虫の脱皮のエクジステロイド代謝に関与するグルコース-メタノール-コリン（GMC）酸化還元酵素のPfamカテゴリーにも機能的に富んでいることがわかった（Supplemental Table S23）。 他の昆虫病原真菌との比較では、O. sinensis系統の遺伝子ファミリーの拡大は、低温適応に関連するとされるPfam用語（PF06772；P < 0.01）の過剰発現でも認められた（Supplemental Table S23）。

一方、収縮状態を示す遺伝子ファミリーは、ABCトランスポーター（PF00005; P < 0.01）, アミノ酸パーミアーゼ（PF00324; P < 0.01）, ATP合成酵素（PF00306; P < 0.05） など輸送過程やエネルギー代謝に関わるものが中心だった（補足 表 S24）。 これらの遺伝子ファミリーの動的進化とは別に、O. sinensisではさらに1,077個（〜13.57％）の種特異的遺伝子が検出された（図2C）。 そのうち318個（〜29.53％）の遺伝子は機能注釈が可能で、デンプン結合（GO：2001070、P < 0.01）、病原性（GO：0009405、P < 0）に関連するGOカテゴリーに有意に濃縮されていた。01）、および細胞壁高分子異化過程（GO：0016998；P < 0.01）（補足表S25）。

真菌病原体の感染を避けるために、昆虫宿主はしばしば活性酸素種（ROS）を急速にたくさん生産して病原体を直接死滅させる。 これに対し、病原体は進化の過程で活性酸素の抗酸化防御システムを進化させ、その中でも活性酸素消去酵素として働くペルオキシダーゼは最も顕著で不可欠な構成要素の一つとみなされている24, 25. O. sinensisで拡張された遺伝子の中で、ペルオキシダーゼ活性は非常に濃縮された機能カテゴリーの一つであった（Supplemental Table S23）。 隠れマルコフモデル（HMM）検索により、O. sinensisのペルオキシダーゼ遺伝子は42個（〜0.53%）であり、この数はC. militaris（28）や酵母（21）よりも著しく多く、ペルオキシダーゼ遺伝子が2倍に拡大することにより、O. sinensisの活性酸素解毒を助ける能力が強くなる可能性があると考えられる（図3Aおよび補足表 S26）。 これら42のペルオキシダーゼ遺伝子のうち、ハロペルオキシダーゼ（haem）が最も多く、検出された全ペルオキシダーゼの約16.67％を占めている（図3B）。 他の近縁真菌種が典型的な2-システインペルオキシレドキシンを完全に欠いているのとは対照的に、O. sinensisはまだ1コピーを保持している（図3B）。 2-Cys peroxiredoxinは、Vibrio vulnificusにおいて、異なるレベルの酸化ストレスに応答する役割を果たすことが以前に示されている26。 比較解析の結果、O. sinensisの保持遺伝子はPrx1型に属し、機能的に保存されており27、外来で生成した高レベルのH2O2に細胞がさらされたときのみ発現すると報告されている26。

植物病原菌（PP）の感染メカニズムは、複雑な植物細胞壁を分解する糖質活性酵素（CAZymes）を必要とするが28、昆虫病原菌（IP）は通常クチクラを貫通し、宿主に感染する29。 このことを検証するために、我々はO. sinensisとさらに4種類の昆虫病原体（M. anisopliae, M. acridum, C. sinensis）を比較した。 militaris、B. bassiana）と4種の植物病原菌（Fusarium graminearum、Magnaporthe grisea、Grosmannia clavigera、Botrytis cinerea）を比較した（補足表S17）。 その結果、昆虫病原菌は植物病原菌に比べて、昆虫のクチクラを分解するプロテアーゼ（IPでは平均362、PPでは342、P < 0.43）およびプロテインキナーゼ（IPでは平均151、PPでは119、P < 0.0014）などを多く有することが示された（補足表 S27-29）。 一方、植物病原菌は昆虫病原菌に比べて、植物細胞壁分解のためのCAZymesを多く保有していた（平均で、IPでは161、PPでは231）（補足表S30〜32）。 植物病原菌を除くと、O. sinensisはM. anisopliae (437), M. acridum (361), C. militaris (355), B. bassiana (396) などの他の昆虫病原菌より、プロテアーゼをコードする遺伝子が著しく少ない（260）ことが分かる。 しかし、O. sinensisで特徴づけられたこれらのプロテアーゼの約35%は、病原体-宿主間の相互作用に関与した可能性が高いシグナルペプチドを含んでおり（補足表S10およびS34）、他の昆虫病原真菌のそれ（平均で20%）よりも大きいことが分かった。 また、他の昆虫病原体と同様に、GH7、GH45、GH51などいくつかのセルラーゼファミリーもO. sinensisでは減少または消失した（補足表S30）。

我々は、菌対虫の長さの比が〜1.20×、〜1.75×、〜2.20×に達するO. sinensisの三つの発生段階にわたる遺伝子発現プロファイルについても検討を行った。 その結果、3つの発生段階間で合計411個の遺伝子が差次的に発現していた（Supplemental Figure S14）。 これらの411個のDEGについて機能アノテーションを行ったところ、グリコシルヒドロラーゼファミリー16（PF00722；FDR < 0.01）、シトクロムP450（PF00067；FDR < 0.01）、メジャーファシリケータースーパーファミリー（PF07690；FDR < 0.05）など主に菌の病原性に関わるものがあることがわかった。 また、NAD dependent epimerase/dehydratase family (PF01370; FDR < 0.01) やBCS1 N-terminal domain (PF08740; FDR < 0.01) などミトコンドリア呼吸鎖に関わる酵素のコード遺伝子も機能的に豊富に含まれていた（補足表 S33）。

Positive Darwinian selection serves as driving forces for fungal pathogenicity

青海チベット高原の高地環境に生息する多様な生物の進化に、正の選択が重要な役割を果たしたことは間違いなく、表現型形質の多くにそのような選択のサインが見られると考えられる3，4，5． O. sinensisと他の12菌種に共通する高信頼性の単一コピーオーソログ1,499個のうち、O. sinensisでは163個の正の選択遺伝子（PSG）が枝サイト尤度比検定（LRT；P < 0.05）により特定できた（補表 S35）。 そのうち1つの遺伝子（OSIN3929；ここではOsPRX1と命名）は、ペルオキシダーゼ活性に機能的に関与していることがわかった（図3C）。 この遺伝子は1-システインを持つペルオキシレドキシンファミリーのメンバーであり、S. cerevisiae 30のPRX1（YBL064C）と高い相同性を有している。 S. cerevisiaeと2つのヒト病原真菌A. fumigatusとC. albicansのPRX1遺伝子は機能的に保存されており、宿主細胞内の酸化バーストの無毒化に必要である31, 32。 特に、イネの病原体として知られる Magnaporthe oryzae の PRX1 を欠損させると、病原性がほぼ完全に消失することから、このペルオキシダーゼが宿主-病原体相互作用の鍵であることが示唆されている27。 また、ペルオキシソーム生合成、プロテインキナーゼ、メタロペプチダーゼなど、宿主-病原体相互作用に関わるいくつかの遺伝子も正の選択を受けていることが検出された（図3C）。

交配システムの進化

子嚢菌では通常交尾型（MAT）遺伝子座により、交尾システムを制御している33。 我々のゲノム解読の結果、O. sinensisはMAT1-2 idiomorphの中にMAT1-2-1の交配型遺伝子を持っているだけでなく、MAT1-1 idiomorphの中に3つの交配型遺伝子（すなわち、MAT1-1-1, MAT1-1-2, MAT1-1-3 ）を持っていた（補図S15B）。 この特徴は、ほぼ全地域にわたる31の自然集団の全ゲノムリシーケンスを用いて検証され、O. sinensisが実際にホモタリンであることが示された（補足図S15Aおよび補足表S36）。 この特徴は、近縁の真菌病原体であるTolypocladium inflatum (MAT1-2)22, C. militaris (MAT1-1)20, B. bassiana (MAT1-1)21, M. anisopliae (MAT1-1)19, M. acridum (MAT1-2)19 が、異型接合型で単一の交配遺伝子のみを持つことと極めて異なるものであった。 よく知られたホモタリン植物病原菌Fusarium graminearum 34と同様に、O. sinensisのこれら2つのMAT遺伝子座の構成から、特にLTRレトロトランスポンに富む特型ゲノム領域内の融合状態が明らかにされた。 ホモタリンO. sinensisとヘテロタリンC. militarisの分裂は174.2MYA近くに起こったと推定され（補足図S13C）、その進化の過程で自己非互換から自己互換への交配システムの変換を何度も受け（補足図S15C）、糸状菌のノイロスポラ属に似ている35。

青海チベット高原の緯度に基づく集団分岐

ゲノム全体の関係と集団分岐を調べるために、青海チベット高原の青海省、四川省、雲南省、甘粛省とチベット自治区を含む既知の分布範囲で、O. sinensisの31アクセッションを収集して再シーケンスしました（補足図S16、補足表S37）。 平均深度10.1×（生データ）のペアエンドリードを合計1億8300万本（配列量約36.68Gb）作成した（Supplemental Table S38）。 これらのデータから、O. sinensisの集団間の関連性を評価するために、816,960個の一塩基多型（SNPs）と48,092個の厳密なインデル（挿入と欠失）のセットを作成しました（補足図S18-19および補足表S39）。 ゲノム変異の大部分（71.1％）は遺伝子間領域にマッピングされ、コーディング領域にマッピングされたサブセット（23.3％、置換率0.86の同義および非同義SNPs 101,997個から成る）があった（補足図S19、補足表S39）。 SNPデータセットを用いて構築された系統樹は、31のアクセッションを低緯度地域から高緯度地域まで、地理的に3つのグループに分けた（図4A）。この結果は、第1および第2固有ベクトルを用いたPCAによっても補強された（図4Bおよび補足図S20A）。 また、祖先集団の推定数（K）を変化させたところ、K=3のとき、3つの異なるグループがPCAと系統樹の再構成と一致した（図4C、補足図S20B）。 低緯度グループのいくつかのアクセッションは、他の2つのグループと比較して、混血の強い証拠を示し、より分散しており、おそらく祖先の多型の共有や最近の導入イベントによる遺伝的多様性が大きいことを示している（図4D,E）。 これらの3つの緯度に基づく集団間の推定集団差分統計量（F ST ）は、低緯度地域、特にチベットのニンギチ地区からの集団の底辺性をさらに明らかにし、それは集団内の塩基多様性（π）が大幅に高く、他の二つの高緯度グループとの集団差分が低いことによってさらに証明された。 (図4C-F）。

さらに、異なるレベルのSNP内容および非同義変異によって影響を受ける遺伝子を調べた（Supplemental Table S40）。 SNP含有量及び／又は非同義変異の多い上位100遺伝子の機能濃縮解析を行ったところ、ポリケチド合成酵素脱水酵素（PF14765；FDR＜3646＞0.01）、KRドメイン（PF08659；FDR＜3646＞0.01）、縮合ドメイン（PF00668；FDR＜3646＞0.01）など主に真菌二次代謝物の代謝に関わることが示された。 また、アシルトランスフェラーゼドメイン（PF00698；FDR < 0.01）、β-ケトアシル合成酵素（PF00109、PF02801；FDR < 0.01）など脂肪酸生合成に関わる機能カテゴリーも充実していた（補足表S41〜42）

．