Comments
Description
Transcript
修士論文 ミドリムシにおける次世代シーケンサを用いた 代謝パスウェイ
NAIST-IS-MT1351022 修士論文 ミドリムシにおける次世代シーケンサを用いた 代謝パスウェイ情報に基づく遺伝子発現解析 奥田 裕樹 2015 年 3 月 4 日 奈良先端科学技術大学院大学 情報科学研究科 情報科学専攻 本論文は奈良先端科学技術大学院大学情報科学研究科に 修士（工学）授与の要件として提出した修士論文である． 奥田 裕樹 審査委員： 金谷 重彦 教授 （主指導教員） 池田 和司 教授 （副指導教官） MD.ALTAF-UL-AMIN 准教授 （副指導教官） 杉浦 忠男 准教授 （副指導教官） 小野 直亮 助教 （副指導教官） ミドリムシにおける次世代シーケンサを用いた 代謝パスウェイ情報に基づく遺伝子発現解析∗ 奥田 裕樹 内容梗概 ミドリムシ Euglena gracilis は嫌気条件下においてバイオディーゼルの燃料とな るワックスエステルを生成することが知られており，石油の代替燃料として期待 されている．本研究では，脂肪酸合成経路を含めた二次代謝経路に関する遺伝子の 同定および環境変動における挙動を調べるために，次世代シーケンサを用いた網 羅的な遺伝子発現解析を行った．はじめに好気条件で培養し，その後嫌気条件下 で 24 時間培養を行ったミドリムシとその葉緑体欠損株の計 8 サンプルにおいて， RNA-Seq および de novo トランスクリプトームアセンブルにより，脂肪酸合成経 路においていくつかの遺伝子発現変動を確認した．また，代謝パスウェイをもとに した発現変動解析によって，協調して変動する遺伝子群を特定した． キーワード 次世代シーケンサ，RNA-Seq，ミドリムシ，代謝パスウェイ ∗ 奈良先端科学技術大学院大学 情報科学研究科 情報科学専攻 修士論文, NAIST-IS-MT1351022, 2015 年 3 月 4 日. i Differential Gene Expression Analysis Based on Metabolic Pathway Information in Engulena gracilis∗ Yuki OKUDA Abstract Euglena gracilis, a unicellular protist that have photosynthesizing chloroplasts produces wax esters under anaerobic condition, many of which can be used as biofuel. RNA-Seq has become one of the main approaches for identifying the differentially expressed genes for various environments and understanding pathways of secondary metabolites in non-model organisms. In this study, transcriptome profiles from neutral and bleached E. gracilis strains with aerobic and anaerobic conditions were investigated using massive parallel sequencing and we detected candidate genes potentially involved in the production of secondary metabolites. We also found coexpressed genes clusters by using pathway-based expression difference analysis. Keywords: Next generation sequencing, RNA-Seq, Euglena gracilis, Metabolic pathway ∗ Master’s Thesis, Department of Information Science, Graduate School of Information Science, Nara Institute of Science and Technology, NAIST-IS-MT1351022, March 4, 2015. ii 目次 目次 iii 図目次 vi 表目次 vii 略語表 viii 第 1 章はじめに 1.1 1 次世代シーケンサによるトランスクリプトーム解析 . . . . . . . . . . . 1 1.1.1 トランスクリプトーム解析を支える技術の変化 . . . . . . . . . 1 1.1.2 RNA-Seq の情報学的側面 . . . . . . . . . . . . . . . . . . . . 2 1.1.3 Reference-based 方式におけるマッピング . . . . . . . . . . . . 3 1.1.4 de novo 方式におけるアセンブル . . . . . . . . . . . . . . . . 5 1.2 ミドリムシについて . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3 本論文の目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.4 本論文の構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 第 2 章 非モデル生物における遺伝子発現解析のための de novo トランスクリ 2.1 プトームアセンブリの性能評価 12 背景 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1.1 非モデル生物における RNA-Seq の戦略 . . . . . . . . . . . . . 12 2.1.2 解析ソフトウェアの問題点 . . . . . . . . . . . . . . . . . . . . 13 2.1.3 サンプルサイズの問題点 . . . . . . . . . . . . . . . . . . . . . 13 iii 2.1.4 ミドリムシにおける RNA-Seq 解析 . . . . . . . . . . . . . . . 14 2.2 本研究の目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 検証に用いるデータセット . . . . . . . . . . . . . . . . . . . . 15 2.3.2 de novo トランスクリプトームアセンブリのソフトウェアの選定 17 2.3.3 評価方法 1：アセンブラごとの性能評価 . . . . . . . . . . . . . 17 2.3.4 評価方法 2：RNA-Seq のデータサイズの評価 . . . . . . . . . . 20 結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 2.5 2.4.1 評価結果 1：アセンブラの結果における統計的指標の評価 . . . 20 2.4.2 評価結果 2：RNA-Seq のデータ量の評価 . . . . . . . . . . . . 23 考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 第 3 章二次代謝経路に関する代謝パスウェイをもとにした発現変動解析 3.1 3.2 背景 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 31 3.1.1 代謝経路と代謝制御 . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.2 代謝パスウェイ . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.3 遺伝子の転写制御 . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.1.4 環境応答による複数遺伝子の転写制御 . . . . . . . . . . . . . . 34 ワックスエステルの起点となる脂肪酸合成経路 . . . . . . . . . . . . . 34 3.2.1 嫌気条件とワックスエステル発酵 . . . . . . . . . . . . . . . . 34 3.2.2 脂肪酸合成経路 . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2.3 desaturase の多様性とミドリムシ特有の高度不飽和脂肪酸合成 経路 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 本研究の目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.4 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.5 3.4.1 RNA-Seq における発現量推定 . . . . . . . . . . . . . . . . . . 39 3.4.2 藻類データベースの作成 . . . . . . . . . . . . . . . . . . . . . 39 3.4.3 パスウェイ解析における発現変動解析 . . . . . . . . . . . . . . 42 結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.5.1 パスウェイ上の発現量推定 . . . . . . . . . . . . . . . . . . . . 45 3.5.2 脂肪酸含有量の変動 . . . . . . . . . . . . . . . . . . . . . . . 50 iv 3.6 考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 第 4 章おわりに 55 謝辞 57 参考文献 58 発表リスト 63 付表 64 v 図目次 1.1 overlap-layout-consensus によるアセンブルアルゴリズムの模式図 7 1.2 de Bruijn グラフによるアセンブルアルゴリズムの模式図 . . . . . 8 1.3 ミドリムシ E. gracilis の顕微鏡写真 . . . . . . . . . . . . . . . . . 10 2.1 ミドリムシのワックスエステル発酵の代謝経路 . . . . . . . . . . . 19 2.2 異なるサンプルサイズにおけるアセンブル結果の統計指標 . . . . . 24 2.3 アセンブルのグラフ構造における分岐とバブルの模式図 . . . . . . 27 2.4 TER における各サンプルサイズのマッピングの可視化 . . . . . . . 29 3.1 生物の持つ desaturase の種類の違い . . . . . . . . . . . . . . . . 3.2 Linoleic acid から Arachidonic acid への代謝経路 . . . . . . . . . 38 3.3 KNApSAcK ファミリーおよび KNApSAcK Bicycle のウェブイ 37 ンターフェイス . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4 遺伝子発現量の MA プロットと p 値のヒストグラム . . . . . . . . 43 3.5 パスウェイ解析における発現変動解析の模式図 . . . . . . . . . . . 44 3.6 ミドリムシにおける短鎖脂肪酸の代謝経路とその発現変動 . . . . . 48 3.7 ミドリムシにおける高度不飽和脂肪酸生合成経路とその発現変動 3.8 遊離脂肪酸の含有量の変化 . . . . . . . . . . . . . . . . . . . . . . 51 3.9 脂肪酸メチルエステルの含有量の変化 . . . . . . . . . . . . . . . . 52 vi . 49 表目次 2.1 サンプルごとのシーケンス結果の統計値 . . . . . . . . . . . . . . . 16 2.2 アセンブリ結果の統計指標 . . . . . . . . . . . . . . . . . . . . . . 21 2.3 アセンブリ結果における既知遺伝子 6 種の相同性検索の対応 . . . . 22 3.1 逆 β 酸化における酵素の発現量変動 . . . . . . . . . . . . . . . . . 46 3.2 アシルキャリアータンパク質に関する酵素の発現量変動 . . . . . . 46 3.3 長鎖脂肪酸に関する酵素の発現量変動 . . . . . . . . . . . . . . . . 47 4.1 藻類データベースに含まれる生物種名と配列数 . . . . . . . . . . . 64 4.1 藻類データベースに含まれる生物種名と配列数 . . . . . . . . . . . 65 4.1 藻類データベースに含まれる生物種名と配列数 . . . . . . . . . . . 66 vii 略語表 RNA-Seq RNA Sequencing bp base pair CoA Coenzyme A MCM methylmalonyl-CoA mutase PNO pyruvate:NADP+ oxidoreductase PDH pyruvate dehydrogenase TER trans-2-enoyl-CoA reductase GAPDH glyceraldehyde-3-phosphate dehydrogenase PEP phosphoenolpyruvate PrDC pyruvate decarboxylase PYK pyruvate kinase KS β-ketosynthase RD β-ketoreductase DH β-hydroxyacyl-dehydrogenase PUFA polyunsaturated fatty acid TCA tricarboxylic acid ATO1 acetyl-CoA thioesterase ECH β-ketoacyl CoA reductase HADH 3-hydroxyacyl-CoA dehydrogenase ECR enoyl-CoA reductase FAS2 fatty acid synthase subunit α FAS1 fatty acid synthase subunit β FABI enoyl-[acp] reductase KAS1 3-oxoacyl-[acp] synthase I FAE1 fatty acid elongase 1 ELO1 elongation of fatty acids protein 1 viii Abbreviation Metabolic Compounds name 和名 C14:0 Myristic acid ミリスチン酸 C16:0 Palmitic acid パルミチン酸 C18:0 Stearic acid ステアリン酸 C18:1 Oleic acid オレイン酸 Linoleic acid リノール酸 C18:3(6,9,12) γ-Linolenic acid γ-リノレン酸 C20:2(11,14) Eicosadienoic acid エイコサジエノイック酸 C20:3(8,11,14) Eicosatrienoic acid エイコサトリエノイック酸 C20:3(8,11,14) Dihomo-γ-linolenic acid ジホモ-γ-リノレン酸 C20:4(8,11,14,17) Eicosatetraenoic acid エイコサテトラエン酸 C20:4(5,8,11,14) Arachidonic acid アラキドン酸 Eicosapentaenoic acid エイコサペンタエン酸 C18:2(9,12) C20:4(5,8,11,14,17) ix 第1章 1.1 1.1.1 はじめに 次世代シーケンサによるトランスクリプトーム解析 トランスクリプトーム解析を支える技術の変化 トランスクリプトームとは，生物が持つ全遺伝子の発現プロファイルのことであ る．遺伝子はタンパク質の設計図であり，タンパク質は細胞内の様々な化学反応を 含めた細胞の機能を規定する物質として，生命活動全般に関わる．すべての生物は 単細胞生物や多細胞生物に限らず生物固有の遺伝子セットを持ちながらも，外部 環境の変化に適応するためや，細胞分化によってその役割を変化させるために，転 写制御機構によって遺伝子の機能する単位を調節している．ゲノム上の遺伝子は， mRNA に転写されタンパク質に翻訳されることよって機能し，これを遺伝子の発 現と呼ぶ．遺伝子の発現は，基本的に遺伝子の周囲に位置する転写制御領域に転写 因子が結合することによって活性または抑制され，mRNA 転写量および翻訳され るタンパク質の総量が変化することによって制御される．細胞中の mRNA の量を 分析することで遺伝子の発現量を推測することを，トランスクリプトーム解析と 呼ぶ． トランスクリプトーム解析はこれまで DNA マイクロアレイを用いた手法が主流 であったが，近年では次世代シーケンサを用いた手法に置き換わりつつある．マ イクロアレイは，チップ上に固定されたプローブと呼ばれる数十 bp の塩基配列が 対応する mRNA と特異的に結合することを利用し，プローブと結合した分子に付 加されたラベルの蛍光を観測することによって，遺伝子の発現量を相対的に測定 する装置である．チップ上に高密度にプローブを配置することによって複数の遺 伝子を並列に測定可能であり，もっとも高密度なものでは数百万種類のプローブ がひとつのチップ上に集積されている．次に，次世代シーケンサを用いたトラン スクリプトームは，細胞内の mRNA の配列を解読し，そこから遺伝子の発現量を 推定する方法である．これを一般に RNA-Seq という．細胞から抽出した mRNA を短い断片にランダムに分割し，その両端にアダプターを付けることでライブラ リを作成する．その後，次世代シーケンサによってアダプター部位から配列の中 央に向かって配列を読み進める．このような配列を読む操作をシーケンスと呼ぶ． 1 次世代シーケンサ（Next generation sequencer）とは，これまで広く用いられて きたサンガーシーケンスとは異なる手法でシーケンスを行う実験機器の総称であ り，第二世代シーケンサ（Second generation sequencer） ，または超並列シーケンサ （Massively parallel sequencer）とも呼ばれる．現在では Illumina 社の HiSeq を 筆頭に，Applied Biosystem 社 (現 Life Technologies) の SOLiD，Roche 社の 454 GS-FLX，Pacific Biosciences 社の PacBio，Life Technologies 社の Ion Torrent などが利用されている． どちらも遺伝子の発現量を計測する手法であるが，マイクロアレイと比較して次 世代シーケンサが持つ利点は，精度の高さと適応範囲の広さにある [1] ．まず精度の 高さとして，次世代シーケンサは測定できる発現レベルのレンジが広いことが挙げ られる．遺伝子発現量はシーケンスされた配列の数によって決定されるため，測定 できる発現量の上限は無く，同時に低発現の遺伝子も高感度で計測が可能である． 一方マイクロアレイは蛍光強度から相対的に発現量を計測するため，次世代シーケ ンサと比較して発現レベルのレンジは狭い [2] ．次に適応範囲の広さとして，マイ クロアレイは原理上，対象とする遺伝子の配列に対応するプローブを作成する必要 がある．そのため，発現量を推定する際には配列が決定されている遺伝子しか扱う ことができず，対象とする生物種も限定される．その点において次世代シーケンサ は，細胞内のトランスクリプトームをシーケンスするため，対象とする生物のゲノ ム配列の情報を必要としない．よって，ゲノム配列が決定されているモデル生物に 加えて，個々の遺伝子の配列や機能の多くが未決定な非モデル生物にも適用するこ とができる．また，発現量推定と同時に，一塩基多型（SNP）や突然変異などの個 体間にまたがる塩基の違いを配列レベルで検出することができるのも大きな特徴で ある． 1.1.2 RNA-Seq の情報学的側面 以上のような背景により，マイクロアレイが利用されてきた研究分野に徐々に次 世代シーケンサが進出しつつある．しかしながら，次世代シーケンサは現状で発展 途中の技術であり，いまだにトランスクリプトーム解析の標準的なワークフローが 確立されていない．そのなかでも特に，蛍光強度を計測値として直接発現量に変換 2 するマイクロアレイとは違い，次世代シーケンサはトランスクリプトームをシーケ ンスすることから，解読された配列データから個々の遺伝子の発現量に変換するた めの計算機的処理が必要である．それは，解読した遺伝子配列がもともとゲノム上 のどの遺伝子から転写されたものなのかを特定する処理である． RNA-Seq における発現量推定の解析は，大きく分けて 2 種類の方法が発案され ている．ゲノム配列をリファレンス配列としてマッピングする Reference-based 方 式と，シーケンスした配列の断片をつなぎ合わせることでリファレンス配列を再構 成する de novo 方式である [3] ． 1.1.3 Reference-based 方式におけるマッピング Reference-based 方式は，シーケンスした断片配列をリファレンス配列にマッピ ングする方法である．トランスクリプトームをシーケンスすることで得られた塩基 配列は，当然ながらその配列が複製された箇所がゲノム上に存在する．その転写が 行われた位置は，配列アラインメントによって対応付けることができる．この操作 をマッピングといい，配列の鋳型として用いるゲノム配列をリファレンス配列と呼 ぶ．シーケンスで得られたすべての配列についてマッピングを行い，対応する断片 配列の数を数えることによって，遺伝子の発現量を推定する．ここで，発現量の値 として代表的なものが，シーケンスで得えられた配列の本数と遺伝子の長さで正 規化した FPKM（Fragments per kilobase of exon per million mapped sequence reads）である．ある遺伝子 t における FPKM は，以下の式によって求められる． FPKM = 109 × Xt lt N ここで，Xt は遺伝子 t にマップされた配列の本数，lt は遺伝子 t の配列長，N は シーケンスで得られた配列の総数である． このような発現量を計算するソフトウェアは，Blat[4] や BWA[5] ，GSNAP[6] ， Bowtie[7] ，Bowtie2[8] など多数公開され，今日に至るまで様々な研究者や企業に よって開発が続けられている．その理由として，情報学的に処理を行うマッピング や発現量推定には幾つもの課題が存在し，ソフトウェアの処理速度などの性能や 発現量推定の結果が大きく異なるからである．そのなかで特に重要なのが，高速 3 にマッピングを行うアルゴリズムと生物が持つ特徴的な転写システムへの対応で ある． 次世代シーケンサは原理的に短い塩基配列を膨大な数シーケンスする．例えば illumina 社の HiSeq 2000 では，およそ 11 日で 600 Gbp ものデータを出力し，マッ ピングに必要な配列数は 60 億本にも及ぶ．それぞれの配列は 100 bp 程度だが，そ れらと一致する箇所をリファレンスとなるゲノム配列の中から探し出さなければい けない．例えばヒトのゲノムサイズは 3.1 Gbp である．ヒトにおける RNA-Seq を 次世代シーケンサ 1 回分実験した場合には，31 億文字の中から文字列が一致する 100 文字を数十億回検索するということになる．このような大規模解析に使われて いるアルゴリズムが，Burrows Wheeler Transform（BWT）である．BWT は別 名ブロックソートとも呼ばれ，もともと文字列の圧縮アルゴリズムとしてよく知ら れている．BWT に加え，ゲノム配列において FM-index と呼ばれる圧縮付全文索 引を作成しておくことにより，効率的かつ高速な部分文字列検索が可能になる． マッピングにおける生物学的な課題として考慮しなければいけないのが，真核生 物におけるスプライシングである．遺伝子は，タンパク質をコードしているエキソ ンと，タンパク質をコードしないイントロンの 2 種類からなり，エキソンの間に複 数のイントロンが挟まる構造をしている．ゲノム上の遺伝子から転写された直後の mRNA 前駆体は，スプライシングという過程においてイントロンが除去され，残っ たエキソンが一つにつなぎ合わさり mRNA となったのちに，翻訳されてタンパク 質となる．そのため，イントロンをまたいだ箇所がシーケンスされた場合，配列と 完全に一致する部分は厳密にはゲノム上に存在せず，2 つの分断された配列として 対応付けられる．また，遺伝子に含まれる複数のエキソンのうち，どれが選ばれる かが場合によって変化するという，選択的スプライシングと呼ばれる現象も起こ る．このとき，ゲノム上の一つの遺伝子配列から，スプライシングバリアントと呼 ばれる数種類の異なる mRNA 配列が転写される．スプライシングされた mRNA では，配列が共通な領域とバリアント同士で配列が異なる領域が混在することにな るため，マッピングによる同定やそれぞれのバリアントの発現量を定量することが 難しくなる．発現量推定を行うソフトウェアは，スプライシングによってマッピン グする配列がリファレンス配列上で分断されることを前提とし，選択的スプライシ ングによって一つの遺伝子の配列から生じる異なった mRNA 由来の配列を認識す 4 る技術が不可欠である．これらの生物学的な課題に対応するには，前述のマッピン グアルゴリズムに加えて，より複雑な部分文字列検索を行う必要がある． 1.1.4 de novo 方式におけるアセンブル de novo 方式は，シーケンサから出力された配列の文字列をもとに，遺伝子配列 を先に再構成する方法である．上述のように次世代シーケンサは遺伝子の断片を大 量にシーケンスするが，配列が読み取られる箇所は完全にランダムである．モデル 生物のように既にリファレンスとなるゲノム配列が解読され遺伝子のアノテーショ ンがされている場合は，次世代シーケンサで読まれた配列を既知のリファレンス配 列と比較することで，その断片配列が何の遺伝子由来かを同定し，その遺伝子の発 現量を定量することができる．一方で，非モデル生物などの参照可能なゲノム配列 がない場合には，まず次世代シーケンサで読まれた配列の断片から元の mRNA 配 列を同定できるようにする必要がある． RNA-Seq の場合，mRNA の発現量が多ければその配列の様々な断片が繰り返 しシーケンスされることになるため，互いに部分的に重複した断片配列を多数得る ことができる．それらの断片の重複した部分を手がかりにつなぎ合わせることによ り，元の遺伝子の配列を復元することが可能である．このような短い配列同士をつ なぎわせて一つの長い配列を作ることをアセンブリと呼び，新規に遺伝子配列の再 構成を行うことを de novo トランスクリプトームアセンブリという． de novo トランスクリプトームアセンブリを実現するアルゴリズムは，ゲノムア センブリで培われた方法をトランスクリプトーム用に改良したものである．代表的 なものに，overlap-layout-consensus と de Bruijn グラフがある [9] ． overlap-layout-consensus によるアセンブル（図 1.1）は，まずシーケンサから得 られた配列をお互いにアラインメントし合い，末端に重複する配列を見つける．そ れぞれの配列を頂点とし，重複する配列を持つ配列同士を辺でつなぐグラフを考え る．その際に，重複する部分配列の文字数を辺の重みとする．そして，グラフの頂 点を 1 回だけ通り，辺の重みの合計が最大になるようになるべく多くの頂点をつ なぐようなパスを探す．そのようにたどった順番に頂点の配列を繋ぎあわせた結果 が，アセンブルした配列として出力される． 5 de Bruijn グラフによるアセンブル（図 1.2）は，シーケンサから得られた配列 を，まず k-mer と呼ばれる k 文字の短い文字列の集合に分割する．まず先頭 1 文 字目から k 文字の文字列，次に 2 文字目から k 文字の文字列，といったようにすべ ての場合で列挙し，k-mer を頂点として隣り合う頂点を辺でつなぐ．これをすべて の配列について行ってグラフを作成し，同じ k-mer の文字を表す頂点を対応させる ことで複数のグラフを結合する．そして，グラフのすべての辺を通るようなパスを 探し，その順番に頂点の配列を繋ぎあわせた結果が，アセンブルした配列となる． overlap-layout-consensus と de Bruijn グラフはどちらもグラフを用いた方法な がらも，両者の一番の違いはグラフからアセンブルした配列への変換効率である． すべての頂点を一度だけ通るようなパスを考える overlap-layout-consensus はハ ミルトンパス問題と呼ばれ NP 困難であるのに対し，すべての辺を通るパスを考 える de Bruijn グラフはオイラーパス問題であり，効率的に解くアルゴリズムが 存在する．そのため，次世代シーケンサが出力するような巨大なデータに対して は，de Bruijn グラフのほうが overlap-layout-consensus よりも計算時間が短く， またアルゴリズムの特性として反復配列などのアセンブルに強い [10] ．その反面， k-mer は A,T,G,C の 4 種類の塩基からなる配列の組み合わせとして，最大 4k 通 りのパターンが考えられるため，k の値を大きくするにつれて必要となる計算機の メモリ量も増大する．de novo トランスクリプトームアセンブリとして開発された Oases[11] ，Rnnotator[12] ，SOAPdenovo-Trans[13] ，Trans-ABySS[14] ，Trinity[15] はすべて de Bruijn グラフを使用している． de novo トランスクリプトームアセンブリにより得られたコンティグの塩基配列 情報から，どのような機能のタンパク質をコードするのかを推定する必要がある． 既知の遺伝子データベースを用いて配列の類似性検索を行うことによって，それぞ れの mRNA 配列が何の遺伝子であるかを対応付ける．このように情報を付与する ことをアノテーションという．また，Reference-based 方式と同様に，アセンブル で得られたコンティグをリファレンス配列としてマッピングを行うことにより，最 終的に遺伝子の発現量が推定できる． 6 overlap-layout-consensus によるアセンブルアルゴリズムの模式図 R1 から R3 の配列が与えられたとき，まずお互いにアラインメントを行う． 図 1.1 重複する配列同士を辺として繋いだグラフを考え，重複文字列数を重みとする． グラフの頂点を一度だけ通り，辺の重みの合計が最大になるようなパスを選択 し，配列をアセンブルする． 7 図 1.2 de Bruijn グラフによるアセンブルアルゴリズムの模式図 R1 から R3 の配列が与えられたとき，まず配列を k-mer に分割したのちに，隣り合う頂 点同士を辺で繋いで de Bruijn グラフを作成する（ここでは k = 3）．すべての 辺を通るようなパスが，アセンブル結果として出力される配列になる． 8 1.2 ミドリムシについて 本研究では，ミドリムシを対象にトランスクリプトーム解析を行った．そこで， ミドリムシについての生物学的特徴を以下に述べる． ミドリムシ Euglena gracilis はユーグレナ目ユーグレナ科ミドリムシ属に属する 微細藻類である．図 1.3 のような外見をしており，植物のように葉緑体を持ち光合 成を行いながらも，鞭毛によって動物のように自発的に運動することができる．ま た，光合成によって栄養を得る独立栄養生物でありながらも，周囲の有機物を取り 込むことで従属栄養生物としても生存可能である． ゲノム配列は解読されていないが，ゲノムサイズはおよそ 3 Gbp と推定されて いる [16] ．これは藻類のモデル生物である Chlamydomonas reinhardtii や Volvox carteri のゲノムサイズの約 30 倍 [17, 18] であるほか，ゲノムが解読されている他の 藻類の 10 から 100 倍程度である．また，ヒトゲノムとほぼ同じ大きさである． ミドリムシは近年様々な応用が期待されている生物である [19] ．現在もっとも有 名な利用方法は食品への添加である．ミドリムシは栄養価が高く，クッキーやラー メンなどに添加され，健康食品として販売されている．また，食用の他にも工業的 利用が期待されており，特定の株のミドリムシは嫌気条件下においてワックスエス テルを合成し体内に蓄積することが知られている [20, 21] ．このワックスエステルは 精製することによってバイオディーゼルの燃料となるため，石油に換わる代替燃料 としての可能性を秘めている． 9 図 1.3 ミドリムシ E. gracilis の顕微鏡写真 “Euglena gracilis” by David Patterson, Linda Amaral Zettler, Mike Peglar and Tom Nerad from Encyclopedia of Life (http://eol.org/data_objects/27470608), is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported (CC BY-NC 3.0). 10 1.3 本論文の目的 そこで本研究では，バイオディーゼルの効率的産生に向けた基礎データを創出す る目的で，ミドリムシの発現プロファイル解析を行う．ミドリムシにおけるワック スエステルの生合成機構を対象に，生育環境の違いによる遺伝子発現の差異を解析 し，嫌気条件下の代謝反応を解析する．また，質量分析器を用いて，好気条件と嫌 気条件での代謝産物の違いを比較する．そして，ワックスエステルの生合成の調節 に関わる代謝経路や遺伝子を見つけ出し，ワックスエステルの生合成を効率よくコ ントロールするための方法を分析することを目的とする． 1.4 本論文の構成 本論文は，次のように構成されている．第 2 章では，de novo トランスクリプ トームアセンブリを行うソフトウェアを選定するための評価を行った．第 3 章では 脂肪酸合成経路に関わる遺伝子の発現差異解析を行うとともに，代謝物質の含有量 を計測することによって，環境変動に対するミドリムシの応答をトランスクリプト ミクスおよびメタボロミクスの両面から解析した．第 4 章では本研究の総括と今後 の課題について述べる． 11 第2章 非モデル生物における遺伝子発現解析のための de novo トランスクリプトームアセンブリの性能評価 2.1 2.1.1 背景 非モデル生物における RNA-Seq の戦略 発現量推定における解析手法の選択は，対象としている生物種のゲノム配列解読 の有無，つまり広義におけるモデル生物か非モデル生物かによって異なる．ヒトや マウスなどのモデル生物の RNA-Seq は，ほぼすべての場合においてゲノム配列を 利用する Reference-based 方式であり，がん細胞のトランスクリプトーム解析など の特殊な場合においてのみ de novo 方式が採用される [22] ．これは前者が精度にお いて圧倒的に高く同時に汎用的だからである．しかしゲノム配列が未決定な非モデ ル生物においては，標準的な選択基準はない．たとえゲノム配列が解読されていた としても，配列の精度やゲノムの長さ，アノテーションなどの付加的な情報が不十 分であり，モデル生物ほどの利点が無い場合が多い． ゲノム配列が未決定な非モデル生物において，Reference-based 方式では近縁種 にゲノム配列が既知の生物が存在する場合に，進化学的に両者の生物が類似した遺 伝配列を持つと仮定する．そして配列のミスマッチやギャップなどのアラインメン トの条件を緩和した上で，マッピングを行う．この方法では，モデル生物のような 近縁種の遺伝子に関する情報を対応付けやすいという利点がある．しかし，進化の 過程において遺伝配列が大きく変わった部分についてはマッピングの精度が下が る．またパラログと呼ばれる進化の過程で生じた遺伝子配列の重複により，同一ゲ ノム内に類似した配列を複数持っている場合には，それらの区別が難しい等の問題 点がある． 一方で de novo 方式は，まずシーケンサから出力された配列情報をもとに mRNA の配列をアセンブルする方法である．対象としている生物が持つ遺伝子配列をア センブルしてリファレンス配列とするために，マッピングの精度が高い．その一方 で，ゲノム配列ではなくトランスクリプトーム配列が持つ独自の特徴として，選択 的スプライシングや遺伝子発現量のばらつきによる配列の再構成の難しさが問題に なる． 12 2.1.2 解析ソフトウェアの問題点 Reference-based 方式と de novo 方式のどちらが優れているかという体系的な評 価方法は現在のところ存在しない [3] ．加えて，それぞれの方式において解析ソフト ウェアが様々な研究機関や企業によって独自に開発され，有償もしくは無償で配布 されている．次世代シーケンサの結果を解析する研究者は，それらの中から適切な ものを選択する必要がある． 2.1.3 サンプルサイズの問題点 また，ソフトウェアの選定以外に，サンプルサイズについても評価が必要になる 場合がある．サンプルサイズとは，1 回の実験において細胞から抽出してシーケン スを行うトランスクリプトームの総量のことである．対象としている生物のトラン スクリプトーム量に対するシーケンス結果のサンプルサイズの比率を，カバレッジ という．カバレッジが 1 倍というのは，対象の生物の持つトランスクリプトームの 全体をちょうど 1 回読み取るだけの量をシーケンスしたということである．カバ レッジが多いほど発現量推定の精度が高くなり，低発現の遺伝子も検出しやすくな る．また，一般的にサンプルサイズが大きくなるほどに，シーケンスに混入するノ イズの量は相対的に減少する．ここで言うノイズとは，シーケンス時の読み取り間 違いであったり（例えば A が T になるなど），配列解読自体の精度が低いことに起 因する読み取り不能な塩基が生じることである（A,T,G,C の代わりに N が挿入さ れる）．そのため，実験計画においてカバレッジを高く設定することが良いとされ るが，1 回の実験でシーケンスできるサンプルの量は決まっている．サンプルサイ ズを増やしてカバレッジを大きくするには，1 回でシーケンスする実験条件の数を 減らして複数回シーケンスを行わなければならない．このように高いカバレッジを 達成しつつ必要な数の実験条件をシーケンスするのは非常に高いコストがかかるた め，カバレッジを必要最低限に抑えて，許容できる範囲内で混入するノイズの比率 を調節する必要がある． Reference-based 方式においては，マッピングにおけるノイズの除去が容易であ る．マッピング時のギャップやミスマッチの許容を厳しくすることで，ノイズの含 まれた配列を意図的に除去することができる．そのために，シーケンスするサンプ 13 ル量が多ければ多いほど精度を上げやすい．一方で，de novo 方式はノイズの除去 が難しい．それはアセンブル結果を評価する方法が存在しないからである．例えば 一塩基の配列の読み取り間違いがあったときに，それがノイズに由来するものなの か，SNP や突然変異による生物個体間の配列多様性に由来するものなのかを区別 することが難しい．そのため，誤った結果を適切にフィルタリングすることができ ない．サンプルサイズを増やしてもノイズを除去できない一方で，サンプルサイズ が少ないとアセンブルが正しく行われない可能性がある．de novo トランスクリプ トームアセンブリにおけるカバレッジとノイズは，トレードオフの関係にある． 2.1.4 ミドリムシにおける RNA-Seq 解析 本研究で対象にしているミドリムシにおいて，次世代シーケンサによる RNA-Seq の報告例はまだ無い．また，ゲノム配列が解読されている系統的に近縁な種も存 在せず，予備実験において Chlamydomonas reinhardtii を対象とした Reference- based 方式によるマッピングでは有効な結果が得られなかった．そのため，ミドリ ムシを対象とした RNA-Seq には de novo 方式による de novo トランスクリプトー ムアセンブリが必要であると判断した．そのために，アセンブリを行うソフトウェ アを選定する必要がある．また，de novo トランスクリプトームアセンブリにおけ る，適切なサンプルサイズを推定してアセンブル結果を評価する必要も生じる． 2.2 本研究の目的 本研究は，次世代シーケンサを用いた RNA-Seq において de novo トランスクリ プトームアセンブリに最適なアセンブラを決定することを目的とする．そのため に，複数のアセンブラを同一条件で実行し，結果の統計指標や遺伝配列の再現度を 解析することで，アセンブリの結果が最も良いアセンブラを検討する．また，サン プルサイズがアセンブリ結果に与える影響を調べるとともに，十分なアセンブルに 必要なサンプルサイズを評価する． 14 2.3 2.3.1 方法 検証に用いるデータセット 本実験には，ミドリムシを対象にした次世代シーケンサによる RNA-Seq のデー タを用いた．ミドリムシの野生株および葉緑体欠損株を対象に，はじめに好気条件 で 2 日または 4 日培養したサンプルと，その後嫌気条件下で 24 時間培養を行った サンプルの合計 8 サンプルについて RNA-Seq を行った．ここで，培養期間 2 日 は対数増殖中期，4 日は対数増殖後期から定常期におけるサンプルである．また， RNA-Seq は Illumina 社の HiSeq 2000 を使用してシーケンスを行った．使用した cDNA ライブラリは 101 bp ずつのペアエンドライブラリであり，調整されたイン サートサイズは 448 bp から 477 bp 程度であった．これらの実験で得られた配列 データの詳細は表 2.1 の通りである． 15 表 2.1 サンプルごとのシーケンス結果の統計値 G は野生株，B は葉緑体欠損 株，2d および 4d は培養期間，an が付与されたものは嫌気条件下で培養した後 のサンプルを表す． Sample ID # of pairs Read length (bp) Total bases (bp) G_2d 20,281,626 101 4.1 × 109 G_2d_an 24,762,047 101 5.0 × 109 G_4d 27,084,024 101 5.5 × 109 G_4d_an 24,194,393 101 4.9 × 109 B_2d 22,594,403 101 4.6 × 109 B_2d_an 24,516,325 101 5.0 × 109 B_4d 19,683,518 101 4.0 × 109 B_4d_an 23,000,158 101 4.6 × 109 16 2.3.2 de novo トランスクリプトームアセンブリのソフトウェアの選定 アセンブリを行うソフトウェアのことを，アセンブラという．今回は無償で利用 できるアセンブラのうち，活発に開発がすすめられており，かつ de novo トランス クリプトームアセンブリによく利用されているアセンブラを選定した．今回の検証 に用いたアセンブラは，Trinity，Oases，SOAPdenovo-Trans の 3 種類である． Trinity[15] は，Broad Institute によって開発された de novo トランスクリプトー ムアセンブラである．このアセンブラは他の 2 つのアセンブラと異なり，トランス クリプトームアセンブリに特化して研究開発が始められた．アセンブリの各工程を Inchworm，Chrysalis，Butterfly という 3 つの独立したソフトウェアに分割し，そ れらが個別に k-mer 頻度をもとにした配列伸長や de Bruijn グラフによるアセンブ ルを順に行っている．本実験では Trinity r2013-02-25 を用いた． Oases[11] は，EMBL European Bioinformatics Institute によって開発された de novo トランスクリプトームアセンブラである．同機関が開発した Velvet と呼ばれ るゲノムアセンブラを基本としており，動作環境に Velvet を必要とする．Velvet によるアセンブルによって基本となるコンティグを生成し，その後に Oases による de Bruijn グラフをもとにしたアセンブルによって最終的な結果を出力する．本実 験では Velvet 1.2.09 および Oases 0.2.08 を用いた． SOAPdenovo-Trans[13] は，Beijing Genomics Institute によって開発された de novo トランスクリプトームアセンブラである．同機関が開発した SOAP および SOAPdenovo と呼ばれるゲノムアセンブラを基本として，Trinity をベースにし た低頻度 k-mer 除去や Oases をベースにしたサブグラフ構築などの機能を盛り 込み，トランスクリプトーム用に改良されたソフトウェアである．本実験では SOAPdenovo-Trans 1.02 を用いた． 2.3.3 評価方法 1：アセンブラごとの性能評価 3 種類の de novo トランスクリプトームアセンブラから出力された結果を，定量 的および定性的に評価した．なお，評価の対象となるアセンブルによって得られた 個々の配列のことを，ここではコンティグと呼ぶ． まず，コンティグの統計指標を算出し比較した．今回用いる統計指標は，コン 17 ティグの総数，コンティグの総塩基長，1,000 bp 以上のコンティグの総数，最大 コンティグ長，N50，平均コンティグ長である．ここで，N50 とはコンティグ長の 加重平均であり，コンティグを長い順に並べて足し合わせた時にコンティグの総 塩基長の半分の長さに到達したときのコンティグ長である．以上の 6 種の指標は， Assemblathon[23] において使用された評価用スクリプト assemblathon_stats.pli を 用いて計算を行った． 次に，既知遺伝子配列との対応付けによって，配列の再現度を評価した．ミドリ ムシにおいて塩基配列が解読されている既知遺伝子の中から，今回の研究において 注目する代謝系を構成する幾つかの遺伝子と，生命維持のために常に一定量発現し ているハウスキーピング遺伝子を選択した．対象とする遺伝子は，ワックスエステ ル発酵の代謝経路に関連する MCM，PNO，PDH および TER をコードする遺伝 子と，ハウスキーピング遺伝子の β-Actin および GAPDH の計 6 種類である．こ れらが機能する代謝経路上の位置は図 2.1 の通りである．また，このハウスキーピ ング遺伝子 2 種は，ミドリムシにおける RT-PCR の発現量正規化に用いられる遺 伝子である [24] ．これらの配列を対象に BLASTx[25] を用いた相同性検索を行い， E-value 1 × 10−6 以下の相同性および配列長の 90% 以上をアセンブルできたかと いう条件で，既知遺伝配列の再現度を評価した． i http://korflab.ucdavis.edu/datasets/Assemblathon/Assemblathon2/Basic_metrics/ assemblathon_stats.pl 18 図 2.1 ミドリムシのワックスエステル発酵の代謝経路 配列の再現度の対象と した遺伝子から生成される酵素を色付きで示した．S. Tucci et al.,2010[21] より 改変． 19 2.3.4 評価方法 2：RNA-Seq のデータサイズの評価 各アセンブラにおいて，サンプルサイズが変化したときのアセンブル結果との関 係を評価方法 1 と同様に評価した．サンプルサイズの評価には，RNA-Seq で得ら れた個別のサンプルの配列データを混ぜあわせた擬似データを用いた．ミドリムシ の野生株から得られた 4 サンプル合計 80 Mbp の配列集合から，非復元抽出による ランダムサンプリングによって 1 Mbp 配列から 80 Mbp 配列までの 10 Mbp 単位 での擬似的な RNA-Seq データを生成した．それぞれのサンプルにおいて，3 種類 のアセンブラでアセンブルをおこなった後に，評価方法 1 と同様のコンティグの統 計指標を算出した． 2.4 2.4.1 結果 評価結果 1：アセンブラの結果における統計的指標の評価 一般にアセンブルの結果の質は，トランスクリプトームアセンブルが網羅的かつ 完全長であることを理想とするため，コンティグの数が多く，それぞれのコンティ グが長いほど良い．3 つの de novo トランスクリプトームアセンブラによって得 られたコンティグの各統計指標を表 2.2 に示した．コンティグの総数やコンティグ の総塩基長，N50，平均コンティグ長の指標において，Trinity がもっとも大きい 値となった．一方で，最大コンティグ長は SOAPdenovo-Trans が最も大きい値と なった． また，アセンブル結果に含まれる既知遺伝子 6 種の有無を調べた結果を表 2.3 に示した．Trinity は全ての既知遺伝子において高い相同性でコンティグが確認 された．Oases はすべての既知遺伝子を確認できたが，MCM において全体の 15%，最長でも GAPDH の 88% と，その一部分しかアセンブルできなかった． SOAPdenovo-Trans は PNO と β-Actin のみにおいて条件を満たすアセンブルが 確認され，他の既知遺伝子も断片的ではあるがすべて 50% 以上の長さのアセンブ ルができていた． 20 表 2.2 アセンブリ結果の統計指標 それぞれの統計指標で最も値が大きい数値 を太字で示している．*は単位が bp を表す． Trinity Oases SOAPdenovo-Trans 210,234 162,632 183,799 170,119,076 65,365,187 66,719,816 # of contigs > 1,000 bp 57,958 9,549 10,970 Longest contig * 12,367 4,752 13,879 N50 * 1,226 550 887 809 402 363 # of contigs Total size of contigs * Mean contig * 21 表 2.3 アセンブリ結果における既知遺伝子 6 種の相同性検索の対応 E-value が小さいほど配列の一致度が高い．また，既知遺伝子の配列長を 100% としたと きのコンティグ長を割合で表している．E-value 1 × 10−6 以下の相同性および 配列長の 90% 以上の両方を満たす数値を太字で示している． Trinity MCM Oases E-value Length(%) 0 100 E-value 1 × 10−72 −138 SOAPdenovo-Trans Length(%) E-value Length(%) 15 0 75 30 0 91 PNO 0 100 1 × 10 PDH 0 100 2 × 10−91 63 0 77 TER 0 95 6 × 10−65 41 0 54 β-Actin 0 100 1 × 10−107 23 0 100 GAPDH 0 100 1 × 10−62 88 1 × 10−151 57 22 2.4.2 評価結果 2：RNA-Seq のデータ量の評価 擬似的にデータ量を変えて同じ実験条件でアセンブルを実行し，サンプルサイズ がアセンブル結果に与える影響を検証した．横軸にサンプルサイズ，縦軸に各統計 指標の値をプロットした結果を図 2.2 に示した．その結果，全体的にコンティグの 総数やコンティグの総塩基長はサンプルサイズの増加にともなって単調に増加し たが，最大コンティグ長や N50，平均コンティグ長においては，サンプルサイズ が 30 Mbp から 40 Mbps 付近で飽和傾向を示した．また，SOAPdenovo-Trans の N50 および平均コンティグ長においては，20 Mbp 付近を境にしてサンプルサイズ の上昇にともない指標の値が減少した． アセンブラによる違いを見ると，最大コンティグ長においてはサンプルサイズ 30 Mbp 以上において SOAPdenovo-Trans が Trinity の結果を上回った．その他 の 5 つの指標においては，Trinity が最も大きい値となった． 23 図 2.2 異なるサンプルサイズにおけるアセンブル結果の統計指標 青の線が Trinity，黄色の線が Oases，赤の線が SOAPdenovo-Trans を示している． 24 2.5 考察 アセンブルされた配列の総数や長さで評価する定量的な尺度においては， Trinity の優位性が確認されたといえる．まず，最大コンティグ長は Trinity と SOAPdenovo-Trans がほぼ同じ結果となったが，1,000 bp 以上のコンティグの総 数および N50 を比較すると，Trinity の結果は全体的に長いコンティグを多く含 む傾向が見られた．一方で，SOAPdenovo-Trans は N50 と平均コンティグの差が 大きいことから，総じて短いコンティグが多く含まれていると考えられる．次に， SOAPdenovo-Trans と Oases は 1,000 bp 以上のコンティグの総数がほぼ同じで あったのにも関わらず，Oases は最大コンティグ長や N50 が低い結果となった．以 上の 2 つの結果から，Oases は他の 2 つのアセンブラと比較して，そもそも長いコ ンティグを生成できなかったとみられる．その理由として考えられるのが，Oases の前処理において Velvet が作成するコンティグにノイズが混入しやすいという ことである．Oases は，Velvet を利用することによって，k-mer から作られる de Bruijn グラフを作るための元の配列を作成する．この段階で配列エラーのノイズが 含まれてしまうことによって，複数の de Bruijn グラフを結合する際にグラフ同士 がうまく繋ぎあわなくなる．その結果として，長いコンティグを生成できなかった 可能性がある．一方で，Trinity と SOAPdenovo-Trans はシーケンサから出力され た配列をそのまま使用し低頻度 k-mer の除去などのノイズ対策を行っているため， 比較的長い配列を作ることができたと考えられる．サンプルサイズを変更した実験 においても，この 2 つのアセンブラはサンプルサイズが 10 Mbp の時点で最長で 10,000 bp 以上の長さの配列がアセンブルできているのに対し，Oases は 70 Mbp でようやく同程度の長さのアセンブルが可能になっている．データサイズが多くな るほど多くのノイズが混入するが，確率的にノイズが全く混入しない正しい配列の 割合も増加し，正しいパスを通るような de Bruijn グラフが作られやすくなる．そ のため，長いコンティグを作るための de Bruijn グラフを作成するためには，多量 のデータが必要であったと推測できる．このことからも，Oases には他の 2 つのア センブラが持つエラー訂正機能の面で顕著な差があると考えられる． 長い配列を多く得られることがアセンブルにおいて良い結果といえる一方で，そ れらの配列がすべて異なる遺伝子由来の配列として生成されるわけではないという ことが考えられる．特に Trinity は，トランスクリプトームの特徴であるスプライ 25 シングバリアントを考慮し，類似した配列を多く出力するという機能がある [26] ． これによって個体間の遺伝子変異や種間においてのスプライシングバリアントの 多様性の解明が期待される一方で，シーケンスされる配列の中にノイズが多く含 まれている場合には，冗長で無意味なコンティグが大量に生成されかねない．つま り，ノイズを含まない正しいコンティグと，ノイズを含む間違ったコンティグの 両方を出力するということである．この問題はアセンブル時のグラフ構造におい て，図 2.3 で示すような分岐（Spur）やバブル（Bubble）が形成されることに起 因し，ヘテロ接合型のゲノムにおけるアセンブルにおいても同様に課題となって いる [9, 27] ．ただし，このような問題の原因となるアセンブラの機能は，Oases や SOAPdenovo-Trans においても同様に存在する．それらがアセンブル結果にどの ように影響を及ぼしているのかは不明であり，スプライスバリアントの存在をゲノ ム配列無しに評価することは難しい．しかしながら，そのようなノイズとなる冗長 性や不要なコンティグは，アセンブル結果の質が良ければアセンブリ以降の解析に よってある程度対処が可能である．具体的には，クラスタリングによって配列冗長 性を除去する方法 [28] や，配列をクラスタリングした上で発現解析を行う方法 [29] ， または単純に既知遺伝子配列の対応付けにおいて対応付けの取れるような質の保証 された配列のみをフィルタリングして使用するという方法がある．これらを実行す るにあたって前提となるアセンブル結果の質を保証するのが，定性的な評価として の既知遺伝子配列との対応付けである．今回の結果から Trinity において特に高い 一致度と完全長で既知遺伝子が検出されたことから，Trinity のアセンブル結果に 占める誤った結果の割合は決して多くないと言える．このように定性的な側面から 定量的評価を支持する結果が得られたことにより，ここで議論した Trinity におけ る配列の冗長性はそれほど問題にならないといえる． 26 図 2.3 アセンブルのグラフ構造における分岐とバブルの模式図 k-mer からな る頂点を四角とし，矢印を辺として重みを矢印の太さで示した．a. 分岐は配列 の末端に配列エラーが生じることによって生じる． b. バブルは配列の中央に配 列エラーが生じることによって生じる．J.R. Miller et al.,2010[9] より改変． 27 また，RNA-Seq のサンプルサイズの評価においては，6 つの評価指標のうちコン ティグ数や総コンティグ長は RNA-Seq のデータ量増加にともなって増加傾向を示 し，最大コンティグ長や N50 などは一定のデータ量で飽和傾向を示すという 2 つ のパターンに大きく分けられた．このような結果は W. R. Francis らのマウスおよ び 6 種の無脊椎動物を対象にした研究 [30] においても同様に確認された．前者の傾 向は，前述のアセンブルにおける冗長性を持つ配列を出力しやすいというアルゴリ ズムによって解釈が可能である．一方で後者の傾向は，特に長い遺伝子が完全長で アセンブルできたことによって，長い配列の割合が増加しなかったと考えられる． サンプルサイズとカバレッジの関係について，図 2.4 は例として評価 1 で使用した TER と各サンプルサイズのマッピングされた配列を可視化している．TER にマッ ピングされた配列の総数とサンプルサイズには強い相関関係があり（r = 0.9998）， 個々の遺伝子配列のカバレッジもサンプルサイズに依存しているといえる．また図 のマッピング結果から，配列長 1,912 bp の TER において配列を先頭から末尾まで を覆うようにシーケンスするのには，少なくとも 10 Mbp から 20 Mbp 以上のサン プルサイズで十分だったことがわかる．配列の取得位置にムラが生じるのは，一般 に塩基配列は GC 含量の違いなどによってシーケンスされやすい部分とシーケンス されにくい部分が存在するからである [31] ．TER は比較的発現が確認されている 遺伝子であるが，他の低発現の遺伝子においてアセンブルを行うには，配列全体を 覆うためにより多くの配列が必要となると思われる．これらの結果を総合すると， おおよそ 20 Mbp から 30 Mbp で十分なカバレッジが得られていると考えられる． 前述の W. R. Francis らの研究においても，20 Mbp から 30 Mbp の範囲において サンプルを作成することがカバレッジとノイズの比率を考慮した際に最良であると いう結果を示しており，本研究の結果と一致する．これらのことから，今回の実験 で取得した各サンプル量の 20 Mbp は妥当であり，アセンブラの選択結果と合わせ て，Trinity において理想的な de novo トランスクリプトームアセンブリが実現で きていると思われる． 28 TER における各サンプルサイズのマッピングの可視化 上から順に 1 Mbp，10 Mbp，20 Mbp，30 Mbp の配列を Bowtie2[8] によりマッピングし IGV[32] で可視化した．40 Mbp 以上のサンプルサイズは，同様の傾向が見られ 図 2.4 たために省略した．灰色の横線が次世代シーケンサから出力された配列を示して いる．灰色の横線が縦に重なるほど，その箇所におけるカバレッジが高いことを 表している． 29 2 つの評価ともに Trinity の優位性を支持する結果が得られたことから，今後の 二次解析以降は Trinity によるアセンブル結果を使用するのが妥当であると判断で きる．また，このようなサンプルサイズの評価結果は，これからの非モデル生物に おけるトランスクリプトーム解析の実験計画を作成する際の基準となりうる．ただ しこれはアセンブル全体を見たときの指標であり，生物学的な要請として，例えば ある発現レベルの遺伝子が最低限取れるためのサンプル量が知りたいといったよう な，具体的な配列取得を保証するものではない．それぞれの実験計画における目標 を達成するには，予備実験をして確認を行うことや実験のレプリケートを多めに設 定するなど，臨機応変な対応が求められる．また，次世代シーケンサの技術的な進 歩はめざましく，Hiseq が出力する一回の実験あたりのシーケンス量は格段に増加 し，PacBio のような平均 10kb 以上のシーケンスが可能なシーケンサも出現してい る．そのような変化に対応するためにアセンブラも改良が繰り返されており，それ らを考慮した総合的な判断が必要である． 30 第3章 二次代謝経路に関する代謝パスウェイをもとに した発現変動解析 3.1 3.1.1 背景 代謝経路と代謝制御 自己複製の基本単位が遺伝ならば，生命活動の基本単位は代謝である．細胞は外 部から取り込んだ物質をもとにして酵素による化学反応を繰り返すことで，エネル ギーを獲得し，生命維持に必要な化合物に変換する．代謝によって生成されたエネ ルギーや化合物は，また別の化学反応に利用される．これら一連の反応の集合は代 謝経路と呼ばれ，物質の変換経路として表される． 生物の恒常性は，代謝経路における代謝物質の流れを適切に調節することによっ て維持される．生体内における代謝は，常に外界の環境変動という形で影響を受け 続ける．周囲の栄養素が枯渇すれば代謝経路の入り口の物質取り込みを行うことが できなくなり，環境中の温度が変化すれば細胞内の酵素活性が変化する．そのよう な外乱に対し，細胞の駆動力としてのエネルギー獲得や，遺伝子の構成単位である 核酸合成系のような一次代謝経路においては，外界の影響を受けにくいロバストな 代謝制御が必要になる．逆に，代謝経路の下流に位置し，個別の機能に特化した二 次代謝経路においては，特定の条件を引き金にして化合物の生成を促進，抑制，遮 断，停止するようなセンシティブな制御が要求される．このような複雑な代謝制御 を実現しているのが，代謝経路の集まりである代謝パスウェイと，化学反応の速度 を規定する酵素である． 3.1.2 代謝パスウェイ 細胞内では数百種類の化学反応により，代謝活動が維持されていると考えられる． それらの化学反応の連鎖を生化学的な機能に基づいてグループ分けし，代謝ネット ワークの部分集合としてまとめたものを代謝パスウェイと呼ぶ．代謝パスウェイの パターンには，直線的な経路もあれば，分岐，ループといった構造を持つ経路もあ る．パスウェイの中には枝分かれする箇所が多数存在し，二つの経路がひとつに合 31 流するなどして，これらが密に組み合わさり連動している．補充経路と呼ばれる迂 回路のようなものも存在し，一次代謝系の定常的な働きに貢献している．代表的な ものに，TCA 回路というループ構造をした代謝パスウェイがある [33] ．これは糖 を分解して得られた acetyl-CoA を消費して，ATP や NADH などの生体内におけ るエネルギー生産に必要な分子を合成する重要な経路である．この回路には 9 個の 中間代謝物が含まれ，化学反応のループを形成している．TCA 回路はエネルギー 生産のための重要な経路として中心的な役割を果たすが，中間代謝物は TCA 回路 以外の代謝の基質としても常に利用される．そのため個々の中間代謝物は，TCA 回路による経路以外からも補充できるようになっており，回路全体のバランスが保 たれる．また，一次代謝においては 1 つの経路が分断されたとしても，他の経路が 働くことによって全体の機能が維持されることがある．大腸菌における研究例で は，解糖系やペントースリン酸回路など関わる遺伝子が欠損している株は，周囲の バイパス経路における酵素量を増やすことで代謝物質量を一定に保っていることが 分かっている [34] ．以上のように，代謝パスウェイは特定の機能を実現するために， 様々なネットワークを構成している．そして，それらはおおまかな機能単位でクラ スタを作成し，互いに綿密に連携することで一つの巨大な代謝経路を形成して生命 活動を支えている． 3.1.3 遺伝子の転写制御 このような複雑な代謝パスウェイの経路を活用して代謝をダイナミックに制御す るのは，反応速度を規定する酵素の物質量，つまり遺伝子の転写量である．酵素は タンパク質によって構成され，ゲノム上の遺伝子から転写および翻訳されることに よって生成される．生体内で起こる生化学反応は，反応の原材料となる基質があれ ば起こるわけではなく，酵素を触媒とすることで初めて反応が始まる．その反応を 調節するために，酵素は重要な役割を果たしている．触媒としての酵素の反応速度 は，酵素自体の濃度や基質濃度，周囲の pH や温度，補酵素や阻害剤などの他のタ ンパク質の存在に大きく依存する．細胞はこれらを巧みに利用することによって酵 素の反応速度を調節し代謝制御を行っているが，基本的な代謝制御を実現するのは 酵素の量である．そして酵素の量は遺伝子の転写によって制御される． 32 転写制御にとってのスイッチのような役割を果たすのが，転写因子である．多く の遺伝子はゲノム上の遺伝配列の前後に制御領域を持ち，その部位に転写因子が結 合することによって，遺伝子発現の促進や抑制が行われる．ゲノム配列の解読がよ り多くの生物種で解明されてきたことによって，より多様な転写因子結合部位とそ の役割が同定されてきた [35] ． 複数の遺伝子を同時に制御するには，一つの転写因子が複数の遺伝子を同じ機構 で発現制御するという構造が最も単純である．このような発現制御機構は，オペロ ンとして古くから知られている [36] ．最初に発見された大腸菌におけるラクトース オペロンは，3 つの異なるタンパク質をコードする領域が隣り合い，その上流に 1 つの制御調節部位を含むプロモーター領域を持つ構造をしている．このオペロンか らは 1 個の mRNA が転写され，別々に翻訳されることによって 3 つのタンパク質 を作る．これらの転写制御には転写因子以外の複数のタンパク質が関与し，制御領 域に結合する制御タンパク質が周囲の他の化合物の濃度によって結合のしやすさを 変化させることで，複雑な発現制御を実現している．このようなオペロンは原核生 物と一部の真核生物で確認されている． オペロンは単純な転写制御の例であり，実際にはより複雑な転写制御が行われて いると考えられる．それは基本的に転写因子に対してタンパク質の量のほうが圧倒 的に多く，転写因子が複数の発現制御に関与しているからである．例えば，植物に おける遺伝子発現解析によって，ストレス耐性を付与する転写因子の 1 つを過剰発 現させた植物体は複数の遺伝子の発現量が変動し，全体としてストレス耐性を上げ るよう代謝を制御していることが示されている [37] ．また，複数の転写因子の組み 合わせによって機能する例として，ヒトの分化後の細胞を iPS 細胞に初期化するた めに，山中ファクターと呼ばれる Oct3/4，Sox2，Klf4，c-Myc の 4 つの転写因子が 協調して働く必要があることは有名である [38] ．複数の転写因子が組み合わさって 機能する例は他にも報告されており [39] ，原理的に転写因子の種類以上の複雑な発 現制御が可能である．出芽酵母においては，3,562 種類の遺伝子と 142 種類の転写 因子の間に 7,074 種類の相互作用があり，転写因子 1 つあたり平均して 50 個の遺 伝子に影響を与えていると推定されている [40] ．このように，遺伝子の発現制御は， 様々な転写因子が協調的に働くことによって，複数の遺伝子の制御が実現される． 33 3.1.4 環境応答による複数遺伝子の転写制御 このように，生物の持つ複雑な代謝制御は代謝パスウェイと遺伝子の転写調節に よって実現しており，互いに切り離せない関係にある．そのような協調関係の具体 例として，代謝パスウェイ上の一連の機能単位において，特定の遺伝子群が協調し て発現することが知られている．これは “guilt-by-association” 原理と呼ばれてい る [41, 42] ．そのなかでも特に二次代謝経路において，パスウェイ上の隣接した複数 の遺伝子が連動して制御されていることが知られている [43, 44, 45] ．このようなパ スウェイ上の遺伝子の発現制御を理解するには，トランスクリプトーム解析による 複数の遺伝子の発現変動を観測するとともに，メタボローム解析による代謝物の定 量が必要になる． 嫌気条件下において，ミドリムシはエネルギー生産の方法を TCA 回路からワッ クスエステル発酵に変更する．この際に，細胞内での発現する遺伝子が大きく変わ ることが知られている． 3.2 3.2.1 ワックスエステルの起点となる脂肪酸合成経路 嫌気条件とワックスエステル発酵 細胞が代謝によってエネルギーを得るには，酸化還元のバランスをとる必要があ る．酸素が十分にある環境下では，酸素を使って糖を酸化して分解し，エネルギー を得ることができる．この過程で酸素分子が還元されて水と二酸化炭素になり，細 胞外に捨てられる．酸素がない環境下では，糖の一部を酸化してエネルギーを得る とともに，残りを還元して細胞外に捨てる．この過程において酢酸を還元すること でエタノールが作られ，β 酸化の経路を逆流することで acetyl-CoA が還元されて 脂肪酸となる．不要となったアルコールは細胞外に捨てられるが，脂肪酸は細胞外 に捨てる代わりにワックスエステルとして細胞内に蓄積される． 34 3.2.2 脂肪酸合成経路 脂肪酸はカルボキシル基を持つ長鎖炭化水素の化合物の総称であり，炭素鎖の長 さと不飽和度の 2 つの組み合わせによって，多様な化合物を構成する．生物が持 つ脂肪酸の多くは直鎖状であり，カルボキシル基から伸びるアルキル基において炭 素鎖が一列に並んだ構造をしている．生体内では acetyl-CoA や malonyl-CoA に よって炭素数が 2 個単位で炭素鎖の伸長が起こるため，炭素鎖が偶数個の偶数脂肪 酸が大半を占める．この炭素鎖の一部は不飽和化されることによって，二重結合に なることがある．二重結合を取ることによって化合物の立体構造に制約が生じ，化 学的な特性が変化する．この二重結合を 2 個以上持つ脂肪酸のことを高度不飽和脂 肪酸（PUFA）という． 脂肪酸は炭素鎖の数と二重結合の数によって示すことができるため，数値とし て表記することがある．例えば，炭素鎖 18 で二重結合を持たない Stearic acid は C18:0，炭素鎖 18 で 1 個の二重結合を持つ Oleic acid は C18:1 と略される． 複数の二重結合を持つ γ-Linoleic acid の場合では，(6,9,12)-Linoleic acid または C18:3(6,9,12) のように二重結合の位置を示す．二重結合の位置は ∆ で表されるた め，(∆6,∆9,∆12)-Linoleic acid のように表記する場合もある． 3.2.3 desaturase の多様性とミドリムシ特有の高度不飽和脂肪酸合成経路 生体内の脂肪酸合成を規定する酵素は，elongase と desaturase の 2 つであ る．elongase は acetyl-CoA や malonyl-CoA を用いて炭素鎖の配列を伸長する． desaturase は炭素鎖の一重結合を二重結合に変化させる．どちらの酵素も炭素 鎖の部位特異的に機能するため，炭素鎖の位置を示すために ∆9 elongase や ∆9 desaturase などと表記される． 脂肪酸はヒトを含め様々な生物の活動を支える重要な化合物であるが，生物 によって desaturase の種類が異なるため，代謝によって合成できる物質が異な る [46] ．例えば，動物は ∆9 以降の炭素を不飽和化できない．代表としてヒト が持つ desaturase は，∆5，∆6，∆9 の 3 種類である．そのため C18:2(9,12) の Eicosapentaenoic acid を始めとして，Linoleic acid や Docosahexaenoic acid の ような高度不飽和脂肪酸を体内で合成できないため，必須脂肪酸として食品から摂 35 取する必要がある．一方で植物や昆虫などは，動物とは違った種類の desaturase を持つ．図 3.1 に示す通り，植物は ∆9，∆12，∆15，昆虫などの無脊椎動物は ∆5， ∆6，∆9，∆12，∆15 の desaturase を持つ． 以上の理由によって，生物はそれぞれに異なる脂肪酸合成経路を持っている．そ の中でも特にミドリムシは，特殊な高度不飽和脂肪酸合成経路を持っている．ミド リムシは多くの生物が持っている ∆6 desaturase が無く，代わりに ∆8 desaturase を持っている．これらは Arachidonic acid や Eicosapentaenoic acid の合成に関与 しているが，代謝経路における反応が ∆6 desaturase から ∆8 desaturase に置き換 わっているわけではなく，elongase による配列伸長と不飽和化の順番が異なってい る．図 3.2 はミドリムシとクラミドモナスにおける Linoleic acid から Arachidonic acid への代謝経路を表しており，お互いに片方の代謝経路しか持っていない． 36 図 3.1 生物の持つ desaturase の種類の違い 脂肪酸の概略図と不飽和化の対 象となる炭素の位置．一般に，動物は ∆5，∆6，∆9，植物は ∆9，∆12，∆15， 昆虫などの無脊椎動物は ∆5，∆6，∆9，∆12，∆15 に特異的に働く desaturase を持つ．Harold W. Cook, 2002[46] より改変． 37 Linoleic acid から Arachidonic acid への代謝経路 クラミドモナ ス Chlamydomonas reinhardtii など一般の生物は，Linoleic acid が不飽和化し た後に伸長反応が起きる．ミドリムシ Euglena gracilis は ∆6 desaturase を持 たず代わりに ∆8 desaturase を持っているため，伸長反応後に不飽和化が起き る．Linoleic acid：C18:2(9,12)，Arachidonic acid：C20:4(5,8,11,14)． 図 3.2 38 3.3 本研究の目的 嫌気条件下におけるミドリムシの環境応答には様々な代謝制御が関わっている と考えられるが，いまだに詳しい機構はよくわかっていない．本章においては， RNA-Seq において得られた遺伝子発現情報とメタボロームの代謝情報を組み合わ せることで，ミドリムシの環境変動に対する応答を明らかにする．そのために，代 謝パスウェイ情報を利用し，複数の遺伝子の協調的な変動を解析する． 3.4 3.4.1 方法 RNA-Seq における発現量推定 ミドリムシにおける RNA-Seq の概要および de novo トランスクリプトームア センブルは 1 章および 2 章に示した通りである．計 8 サンプルの RNA-Seq デー タから Trinity を用いて de novo トランスクリプトームアセンブルを行い，その結 果 208,207 のコンティグを得た．その後 Usearch clustering[28] により冗長な配列 を除去し，122,864 の配列を得た．この配列をリファレンス配列とし，Bowtie[7] を 用いてマッピングを行い FPKM として発現量を推定した．各サンプルの発現値は サンプルサイズのバイアスを除去するために，quantile 正規化法を用いて正規化し た．その結果，83,143 の配列の発現値が得られた．また，これらの配列は BLASTx を用いて KNApSAcK データベースおよび KEGG データベースとの対応付けを行 い，アノテーション情報を付与した． 3.4.2 藻類データベースの作成 RNA-Seq のアセンブルで得られた配列は，独自に作成した藻類データベースに よってアノテーション付けを行った．このデータベースは KNApSAcK Bicyclei と して公開した． KNApSAcK Bicycle は，藻類と分類される生物種のみを対象にしたデータベース である．NCBI の研究者が手動でクオリティをチェックし維持・管理しているデー i http://kanaya.naist.jp/bicycle/top.html 39 タベースである RefSeq[47] の中から，藻類と分類される生物種に該当する計 96,641 のアミノ酸配列とそのメタデータを収集した．そして，各配列に EC Number と Reaction を手動でアノテーション付けした．なお，対象とした生物群はハプト 藻（Haptophyceae） ，クリプト藻（Cryptophyta） ，渦鞭毛藻（Stramenopiles） ，不 等毛藻（Dinophyceae），緑藻（Green algae），紅藻（Red algae），ユーグレナ藻 （Euglenaceae）であり，対象とした生物種の学名と取得した配列数は付表の通りで ある． KNApSAcK Bicycle データベースは，他の KNApSAcK ファミリーと同様に， ウェブインターフェイスからキーワードや配列相同性による検索が可能である． KNApSAcK ファミリーのウェブインターフェイスと，KNApSAcK Bicycle にお ける検索例を図 3.3 に示した．KNApSAcK Bicycle の検索結果において，特定の 対応する化合物についてはリンクが貼れれており，他の KNApSAcK ファミリーの データベースの中から対応する化合物を参照することができる． 40 図 3.3 KNApSAcK ファミリーおよび KNApSAcK Bicycle のウェブイ ンターフェイス a. KNApSAcK ファミリーの表示画面．画面右下に KNAp- SAcK Bicycle へのリンクがある．b. KNApSAcK Bicycle の検索結果の画 面．ここでは AGID:AG0082104 の NADH:ubiquinone oxidoreductase subunit 4L の結果を表示している．Equation カラムにおいて，この化合物が反応 を触媒する基質である NADH と ubiquinone にリンクが貼れらており，他の KNApSAcK ファミリーの中から対応するエントリを参照することができる． 41 3.4.3 パスウェイ解析における発現変動解析 細胞の代謝経路において化学反応が効率よく進むには，一連の化学反応が連鎖し て起こる必要がある．そのため，関連する酵素は共同して発現制御を受けている場 合が多い．すなわち，反応経路の起点となる反応を触媒する酵素が発現していると きには，その下流の酵素も同時に発現している可能性が高い． 個々の遺伝子において t 検定を用いた発現変動解析を行った結果，False Discov- ery Rate（FDR）が 0.05 以下の有意に発現している遺伝子は合計 3,433 個であっ た．図 3.4 に各遺伝子の嫌気条件前後での発現変動および p 値のヒストグラムを 示す．このように，個々の遺伝子においてはそれぞれの発現変動が検出されたもの の，脂肪酸合成経路などにおける代謝パスウェイ上の遺伝子が同時に発現を変動さ せている結果は得られなかった．DNA マイクロアレイや RNA-Seq における遺伝 子の発現量の測定は，どうしても実験による測定誤差が含まれる．今回のようなサ ンプル数の少ない実験においては，遺伝子の発現量の変動を正規化することができ ない場合がある．異なる条件下で遺伝子の発現変動を比較する場合には，機能的に 関連する遺伝子をセットとして統計検定を行うという手法が広く使われている． そこで，協調して発現している遺伝子を推定するために，代謝パスウェイ上で隣 接している任意の n 個の遺伝子のセットについて発現変動解析を行う．本研究で使 用する代謝パスウェイを考慮した発現変動解析の模式図を図 3.5 に示す．例えば， n = 2 ならば，EZM1 と EZM2 のように，隣り合う 2 つの遺伝子を一組にして発 現が変動しているかを検定する．これをすべての組み合わせにおいて行う．また， n = 3 と n = 4 においても同様に検定する． 本研究では野生株と葉緑体欠損株について 4 つのサンプルが存在するため，同 じ実験条件における違う株について遺伝子発現が有意に変動しているかをマン・ ホイットニーの U 検定で検定した．データ処理および統計検定は GNU Rii を使用 した． ii http://cran.r-project.org/ 42 図 3.4 遺伝子発現量の MA プロットと p 値のヒストグラム a. 遺伝子発現変 動の MA プロット．ここでは例として葉緑体欠損株における嫌気条件前後の 2 サンプルずつの比較を示している．各遺伝子の発現量において，横軸が log 2 変 換した発現量の平均（A），縦軸が log 2 変換された発現量の差（M）を表す．A が大きいほど発現量が高く，M が正に大きいほど嫌気条件前で高発現，負に大き いほど嫌気条件後で高発現となる．有意に発現変動がみられた遺伝子を赤で強調 した． b. p 値のヒストグラム．FDR が 0.05 以下である遺伝子を有意な発現変 動がみられた遺伝子として抽出した． 43 図 3.5 パスウェイ解析における発現変動解析の模式図 CPD は代謝物質， EZM は酵素（遺伝子）を表す．隣り合う n 個の遺伝子の組み合わせについてそ れぞれ検定を行う． 44 3.5 3.5.1 結果 パスウェイ上の発現量推定 ワックスエステル合成経路における嫌気条件前後の酵素の発現量を表 3.1，表 3.2，表 3.3 に示す．また，遺伝子の発現量変動を付与した代謝パスウェイについ て，図 3.6 を短鎖脂肪酸の合成経路，図 3.7 を高度不飽和脂肪酸生合成経路として 示した． まず，短鎖脂肪酸の合成経路は図 3.6 のようになっており，上部に示す acyl-CoA を起点として，左側の経路がミドリムシに特異的な逆 β 酸化，右側の経路が一般的 に用いられるアシルキャリアータンパク質を用いた経路となっている．逆 β 酸化に 関わる ATO1，ECH/HADH，ECR はともに有意に発現上昇がみられたが，アシル キャリアータンパク質を用いた経路においては有意な発現変動が見られなかった． 次に，高度不飽和脂肪酸の合成経路は図 3.7 のようになっており，炭素鎖 14 で 不飽和度 0 の Myristic acid（C14:0）を起点として，炭素の伸長および不飽和化に よって最終的に Eicosapentaenoic acid（C20:5(5,8,11,14,17)）が合成される．この 生合成経路では，後半に関与する ∆5，∆8，∆12 desaturase および ∆9 elongase に おいて有意な発現上昇がみられた．その反面，生合成経路の前半に関与する KAS1， FAE1，∆9 desaturase において有意に発現減少がみられた． 45 表 3.1 逆 β 酸化における酵素の発現量変動 それぞれの発現量は同条件の 4 サ ンプルの平均を表す．0h は嫌気条件下で培養する前のサンプル，24h は嫌気条 件下で 24 時間培養後のサンプルを表す． CoA based fatty acid synthesis 0h 24 h acetyl-CoA acyltransferase 9.0 11.3 β-ketoacyl CoA reductase / 3-hydroxyacyl CoA dehydratase 14.6 15.7 enoyl-CoA reductase 10.3 14.0 表 3.2 アシルキャリアータンパク質に関する酵素の発現量変動 それぞれの発 現量は同条件の 4 サンプルの平均を表す．0h は嫌気条件下で培養する前のサン プル，24h は嫌気条件下で 24 時間培養後のサンプルを表す． Acyl carrier protein(acp) based fatty acid synthesis 0h 24 h β-ketoacyl-[acp] synthase 57.8 49.3 β-ketoacyl-[acp] reductase 44.2 43.5 3-hydroxyacyl-[acp] dehydratase 34.2 26.3 enoyl-[acp] reductase 12.8 17.3 46 表 3.3 長鎖脂肪酸に関する酵素の発現量変動 それぞれの発現量は同条件の 4 サンプルの平均を表す．0h は嫌気条件下で培養する前のサンプル，24h は嫌気 条件下で 24 時間培養後のサンプルを表す． Long chain fatty acid synthesis 0h 24 h KAS1 57.8 49.3 FAE1 40.3 39.4 ELO1 18.6 19.7 ∆9 desaturase 41.6 24.5 ∆12 desaturase 150.8 164.6 ∆9 elongase 151.8 140.8 ∆8 desaturase 80.9 93.5 ∆5 desaturase 18.0 23.0 ∆17 desaturase 73.5 67.5 47 図 3.6 ミドリムシにおける短鎖脂肪酸の代謝経路とその発現変動 左の経路が ミドリムシに特徴的な系，右の経路が生物一般に広く用いられる系である．発現 上昇が確認された酵素名を青，発現変動が見られなかった酵素名を灰色で示した． 48 図 3.7 ミドリムシにおける高度不飽和脂肪酸生合成経路とその発現変動 発現 上昇が確認された酵素名を青，発現減少が確認された酵素名を赤，発現変動が見 られなかった酵素名を灰色で示した． 49 3.5.2 脂肪酸含有量の変動 トランスクリプトーム解析によって検出された酵素の発現変動の影響を見るため に，これらが反応を触媒する脂肪酸の含有量を嫌気条件前後で測定した．図 3.8 お よび図 3.9 は，遊離脂肪酸および脂肪酸メチルエステルの含有量を示している．遊 離脂肪酸およびメチルエステルの含有量は，ともに総量に対する相対量である．ガ スクロマトグラフ質量分析計（GC-MS）での測定結果から，C10:0 から C24 まで それぞれのピークの面積値を求め，その合計で割った値をプロットしている．野生 株および葉緑体欠損株において，Myristic acid（C14:0）の含有量は有意に上昇し ており，Palmitic acid（C16:0）と Stearic acid（C18:0）の含有量は減少していた． 50 図 3.8 遊離脂肪酸の含有量の変化 嫌気条件前と後 (0h および 24h) における 各脂肪酸の含有量を示した．上図が野生株，下図が葉緑体欠損株のものである． エラーバーは 3 回の実験結果の標準偏差を示す． 51 図 3.9 脂肪酸メチルエステルの含有量の変化 嫌気条件前と後 (0h および 24h) における各脂肪酸の含有量を示した．上図が野生株，下図が葉緑体欠損株のもの である．エラーバーは 3 回の実験結果の標準偏差を示す． 52 3.6 考察 以上のトランスクリプトーム解析およびメタボローム解析の結果を総合する．嫌 気条件下で 24 時間培養したミドリムシにおいて，Myristic acid が含有量が増加し た．これは短鎖脂肪酸の合成経路における ATO1，ECH/HADH，ECR の発現量 が増加したことによって，より多く合成が進んだと考えられる．また，Stearic acid および Palmitic acid は含有量が減少した．これは高度不飽和脂肪酸の生合成経路 前半に関わる KAS1，FAE1，∆9 desaturase の発現量が減少したことによって合 成速度が減少し，同時に Stearic acid 以降の ∆5，∆8，∆12 desaturase および ∆9 elongase の発現量が上昇したことによって，下流の反応速度が上昇し代謝が促進さ れたためであると考えられる．これらの結果から，ミドリムシは嫌気条件下におい て炭素鎖 14 の Myristic acid を蓄積し，炭素鎖 18 以降の各種脂肪酸の合成を抑制 する傾向にあることがわかった． 本研究が対象としているバイオディーゼルは，炭素数の異なる脂肪酸の種類や配 分量の違いにより，燃料の用途に合わせて揮発性などの性質を調節している．バイ オディーゼル産生の効率化には，それぞれの目的に合わせて代謝経路情報をもとに した必要な化合物産生を促進する必要がある．例えば炭素数の多い脂肪酸を効率よ く合成するためには，Myristic acid 以降の代謝経路がボトルネックとなると考え られる．このような問題に対処するために，特定の代謝経路の遺伝子を過剰発現さ せてボトルネックとなる代謝反応を意図的に促進するような遺伝子改変を行うこと で，炭素数の多い脂肪酸から精製されるバイオディーゼル生産の効率化が実現でき ると考えられる．また炭素数の少ない脂肪酸においても，本研究から代謝経路に関 与する遺伝子が全体的に発現上昇していることが確認された．これらの遺伝子をよ り機能させるような遺伝子改変を行うことで，同様に炭素数の少ない脂肪酸から精 製されるバイオディーゼルの効率的産生が可能になると考えられる． 嫌気条件の有無は脂肪酸合成経路に影響を及ぼした一方で，野生株と葉緑体欠損 株，培養期間の差においては，大きな変化が見られなかった．以下にその理由を考 察する． 野生株と葉緑体欠損株については，そもそも光合成の能力の有無が脂肪酸合成と 直接は関係が無いことが考えられるほか，ミドリムシの栄養貯蔵が影響している可 能性がある．植物は一般に光合成によってデンプンを根などに蓄積するが，ミドリ 53 ムシは光合成によって細胞内にパラミロンと呼ばれる多糖類を蓄積する．葉緑体欠 損株においても，従属栄養生物として外界から栄養を摂取してパラミロンを生成す る．この栄養源をワックスエステル発酵に利用することで，野生株と葉緑体欠損株 のどちらでもエネルギー生産が可能になる．ミドリムシのワックスエステル発酵が 外部の二酸化炭素を炭素源としているのか，体内のパラミロンを炭素源にしている のかは不明である．生体内の炭素源の流れを調べるには，代謝パスウェイの解析で は不十分であり，13 C 安定同位体を標識とした動的代謝プロファイリングが必要と なる．少なくとも今回の実験においては，野生株と葉緑体欠損株のワックスエステ ルに至る経路の差を検出できなかった． また，培養期間の差は，対数増殖中期もしくは対数増殖後期から定常期を表す． 増殖期の差が遺伝子発現や代謝物質の含有量に影響しなかったことから，増殖期の 差は嫌気条件の差よりも軽微であることがわかる．ワックスエステル発酵による嫌 気条件下のエネルギー獲得は，一般的な好気条件下のエネルギー生産よりも効率が 悪いため，増殖にエネルギーを利用できなかったと考えられる．しかしながら，増 殖速度は環境中の細胞密度や栄養源の多寡など様々な要因に影響されるため，実際 の増殖の有無を培養時間から推測することは困難である．実際に細胞数の計測など を行うなどの比較をしなければ，増殖期の差というものを明確に捉えられない．結 果的に，今回の場合では 2 日と 4 日という培養期間を大きく変えていて，なおかつ 培養期間の差は小さかったことから，増殖という活動自体はワックスエステルの生 合成にあまり影響を及ぼさないといえる． 今回の研究により，ワックスエステルの生合成に至る代謝パスウェイにおいて， 協調して発現する遺伝子を同定することができた．また，同時にメタボローム解析 の情報を用いることで，発現量変動が及ぼす影響を含めて総合的に解釈することが 可能となった．以上のように，代謝パスウェイを利用した遺伝子の発現変動解析は， 環境変動に応答する特定の代謝経路を解析するのに有用である．また，本研究の手 法は，エンリッチメント解析（Gene Set Enrichment Analysis）などの遺伝子オン トロジー情報を用いた解析よりも，少ないサンプル数で実行することができる．ま たこの手法は，遺伝子オントロジー情報が不十分な非モデル生物においても利用可 能と考えられる． 54 第4章 おわりに 本論文で，次世代シーケンサにおけるトランスクリプトーム解析により，ミドリ ムシの環境変動に対する遺伝子発現変動を解析した． 第 2 章では，de novo トランスクリプトームアセンブリにおけるソフトウェアの 性能を評価し，アセンブリに最適なサンプルサイズを検討した．選出した 3 種類の de novo トランスクリプトームアセンブラの中では Trinity が最も優れた結果とな り，サンプルサイズの多寡に関わらず安定した結果を示した．Google Scholari にお ける 2015 年 2 月 1 日現在の Trinity の論文引用数は 1,410 と，Oases の 372 およ び SOAPdenovo-Trans の 35 を大きく引き離している．このことからも，Trinity が de novo トランスクリプトームアセンブラとして研究者に支持されているといえ る．本研究で Tirnity の性能の優位性を定量的および定性的に評価したことは，今 後の de novo トランスクリプトームアセンブリを用いた研究におけるアセンブラ選 択の根拠として貢献できると考える． 第 3 章では，トランスクリプトーム解析とメタボローム解析の両面から，環境応 答に対する代謝と発現制御の関係性について解析した．その結果，ワックスエステ ル発酵からワックスエステル合成に至るまでの経路において，遺伝子の制御が行わ れ，代謝物質の含有量が変化していることがわかった．本研究では，直接ワックス エステルの収量増加といった応用分野には踏み込めていない．しかしながら，代謝 経路を同定して代謝物質の含量の変動を示すことができたほか，関連する遺伝子の 配列を取得することができたことから，基礎研究として成果を上げられたと考えら れる．また，本研究で対象としたミドリムシのワックスエステルを利用した代替燃 料の生成以外にも，大腸菌やクロストリジウム属にミドリムシ由来の TER を発現 させることで，工業利用可能な 1-butanol の生産効率を上げることに成功している 事例がある [48] ．ミドリムシの特殊な遺伝子や代謝経路が持つ潜在的な利用価値は 高い．ミドリムシの環境応答に関する知見は，これからのミドリムシの工業的利用 や有用な合成経路の同定の研究に貢献できると考えられる． トランスクリプトーム解析は，次世代シーケンサが標準的な実験手法となる近い 将来においても，非モデル生物を含めたマイノリティの研究対象において大きく力 i https://scholar.google.com 55 を発揮する手法であると考えられる．トランスクリプトーム解析において，対象と する生物のゲノム配列が利用できるということは，大いに研究の幅が広がることを 意味する．しかし，次世代シーケンサ登場以降，様々な種においてゲノム配列の決 定が相次いで報告されているが，いまだにドラフトゲノムとして未完成なものが多 い．モデル生物のゲノム配列ですら，今なお完全な解読に向けて情報が更新され続 けていることを考えると，非モデル生物においてモデル生物が辿ってきた研究の歴 史を追従する場合，莫大な時間とコストがかかることは明白である．次世代シーケ ンサを用いた研究におけるトランスクリプトーム解析は，ゲノム解読の代替案とし て非常に有効な選択肢である．これからのデータドリブンな生物学において，トラ ンスクリプトームのデータが蓄積されていくことにより，ゲノム配列中心の研究ス タイルが一変する可能性がある．そのような分野において，情報学的な側面から次 世代シーケンサの解析手法を研究することは大きな意味を持つと考えられる． 56 謝辞 主指導教官の金谷 重彦先生には，本研究室で研究する機会を与えてくださったこ とに感謝致します．研究者としてのアイデアと想像力の大切さ，何事にも怖じけず 試してみようとするチャレンジ精神，そして活力こそ研究におけるすべての源泉な のだという姿勢を学ぶことができました．ここに感謝致します． 副指導教員の池田 和司先生には，ゼミナール II において貴重な助言を頂きまし た．ここに感謝致します． 副指導教員の MD.ALTAF-UL-AMIN 先生および杉浦 忠男先生，佐藤 哲大先 生，西岡 孝明先生，黄 銘先生には，研究室のゼミナールにおいて多くの助言をい ただき，また指導してくださいました．ここに感謝致します． 副指導教員の小野 直亮先生には，本研究の指導を始めとして，研究室の活動や自 主的な勉強会などを含め，様々な活動においてご助力頂きました．ここに感謝致し ます． 国立遺伝学研究所 生命情報研究センター 有田 正規先生および大阪府立大学 生 命環境科学研究科 太田 大策先生には，本研究を進めるにあたって多くのご援助を 頂きました．ここに感謝致します． 本学の在籍期間中に支えてくださった計算システムズ生物学研究室の皆様に感謝 致します． そして最後に，私を情報学とバイオインフォマティクスという研究分野に導いて くださった，ライフサイエンス統合データベースセンター 坊農 秀雅氏， 産業技術 総合研究所 ゲノム情報研究センター 瀬々 潤氏に感謝致します．本当にありがとう ございました． 57 参考文献 [1] Zhong Wang, Mark Gerstein, and Michael Snyder. “RNA-Seq: a revolutionary tool for transcriptomics.” Nature Reviews Genetics 10.1 (2009): 57-63. [2] Ali Mortazavi, et al. “Mapping and quantifying mammalian transcriptomes by RNA-Seq.” Nature methods 5.7 (2008): 621-628. [3] Jeffrey A. Martin and Zhong Wang. “Next-generation transcriptome assembly.” Nature Reviews Genetics 12.10 (2011): 671-682. [4] W. James Kent. “BLAT―the BLAST-like alignment tool.” Genome research 12.4 (2002): 656-664. [5] Heng Li and Richard Durbin. “Fast and accurate short read alignment with Burrows–Wheeler transform.” Bioinformatics 25.14 (2009): 1754-1760. [6] Thomas D. Wu and Serban Nacu. “Fast and SNP-tolerant detection of complex variants and splicing in short reads.” Bioinformatics 26.7 (2010): 873-881. [7] Ben Langmead, et al. “Ultrafast and memory-efficient alignment of short DNA sequences to the human genome.” Genome Biology 10.3 (2009): R25. [8] Ben Langmead and Steven L. Salzberg. “Fast gapped-read alignment with Bowtie 2.” Nature methods 9.4 (2012): 357-359. [9] Jason R. Miller, Sergey Koren, and Granger Sutton. “Assembly algorithms for next-generation sequencing data.” Genomics 95.6 (2010): 315-327. [10] Pavel A. Pevzner, Haixu Tang, and Michael S. Waterman. “An Eulerian path approach to DNA fragment assembly.” Proceedings of the National Academy of Sciences 98.17 (2001): 9748-9753. [11] Marcel H. Schulz, et al. “Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels.” Bioinformatics 28.8 (2012): 10861092. [12] Jeffrey Martin, et al. “Rnnotator: an automated de novo transcriptome assembly pipeline from stranded RNA-Seq reads.” BMC genomics 11.1 (2010): 663. 58 [13] Yinlong Xie, et al. “SOAPdenovo-Trans: de novo transcriptome assembly with short RNA-Seq reads.” Bioinformatics 30.12 (2014): 1660-1666. [14] Gordon Robertson, et al. “de novo assembly and analysis of RNA-seq data.” Nature methods 7.11 (2010): 909-912. [15] Manfred G. Grabherr, et al. “Full-length transcriptome assembly from RNA-Seq data without a reference genome.” Nature biotechnology 29.7 (2011): 644-652. [16] Muh-Ching Yee and James C. Bartholomew. “Light regulation of the cell cycle in Euglena gracilis bacillaris.” Cytometry 9.4 (1988): 387-393. [17] Sabeeha S. Merchant, et al. “The Chlamydomonas genome reveals the evolution of key animal and plant functions.” Science 318.5848 (2007): 245-250. [18] Simon E. Prochnik, et al. “Genomic analysis of organismal complexity in the multicellular green alga Volvox carteri.” Science 329.5988 (2010): 223226. [19] 亮 嵐田. “微細藻類ユーグレナの特徴と食品・環境分野への応用.” Journal of the Society of Tropical Resources Technologists 26.1 (2010): 19-22. [20] Hiroshi Inui, et al. “Wax ester fermentation in Euglena gracilis.” FEBS Letters 150.1 (1982): 89-93. [21] Sara Tucci, et al. “Variability of wax ester fermentation in natural and bleached Euglena gracilis strains in response to oxygen and the elongase inhibitor flufenacet.” Journal of Eukaryotic Microbiology 57.1 (2010): 6369. [22] Jeyanthy Eswaran, et al. “RNA sequencing of cancer reveals novel splicing alterations.” Scientific reports 3 (2013). [23] Dent Earl, et al. “Assemblathon 1: a competitive assessment of de novo short read assembly methods.” Genome research 21.12 (2011): 2224-2241. [24] Claudine Deloménie, et al. “A new homolog of FocA transporters identified in cadmium-resistant Euglena gracilis.” Biochemical and biophysical research communications 358.2 (2007): 455-461. [25] Stephen F. Altschul, et al. “Basic local alignment search tool.” Journal of 59 molecular biology 215.3 (1990): 403-410. [26] Brian J. Haas, et al. “de novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis.” Nature protocols 8.8 (2013): 1494-1512. [27] Rei Kajitani, et al. “Efficient de novo assembly of highly heterozygous genomes from whole-genome shotgun short reads.” Genome research 24.8 (2014): 1384-1395. [28] Robert C. Edgar. “Search and clustering orders of magnitude faster than BLAST.” Bioinformatics 26.19 (2010): 2460-2461. [29] Nadia M. Davidson and Alicia Oshlack. “Corset: enabling differential gene expression analysis for de novo assembled transcriptomes.” Genome biology 15.7 (2014): 410. [30] Warren R. Francis, et al. “A comparison across non-model animals suggests an optimal sequencing depth for de novo transcriptome assembly.” BMC genomics 14.1 (2013): 167. [31] Yuval Benjamini and Terence P. Speed. “Summarizing and correcting the GC content bias in high-throughput sequencing.” Nucleic acids research (2012): gks001. [32] James T. Robinson, et al. “Integrative genomics viewer.” Nature biotechnology 29.1 (2011): 24-26. [33] Oliver E. Owen, Satish C. Kalhan, and Richard W. Hanson. “The key role of anaplerosis and cataplerosis for citric acid cycle function.” Journal of Biological Chemistry 277.34 (2002): 30409-30412. [34] Nobuyoshi Ishii, et al. “Multiple high-throughput analyses monitor the response of E. coli to perturbations.” Science 316.5824 (2007): 593-597. [35] Michael Levine and Robert Tjian. “Transcription regulation and animal diversity.” Nature 424.6945 (2003): 147-151. [36] François Jacob, et al. “Operon: a group of genes with the expression coordinated by an operator.” CR Acad. Sci. 250 (1960): 1727–1729. [37] Kyonoshin Maruyama, et al. “Metabolic pathways involved in cold accli- 60 mation identified by integrated analysis of metabolites and transcripts regulated by DREB1A and DREB2A.” Plant physiology 150.4 (2009): 19721980. [38] Kazutoshi Takahashi and Shinya Yamanaka. “Induction of pluripotent stem cells from mouse embryonic and adult fibroblast cultures by defined factors.” Cell 126.4 (2006): 663-676. [39] Aika Terada, et al. “Statistical significance of combinatorial regulations.” Proceedings of the National Academy of Sciences 110.32 (2013): 1299613001. [40] Nicholas M. Luscombe, et al. “Genomic analysis of regulatory network dynamics reveals large topological changes.” Nature 431.7006 (2004): 308312. [41] Kazuki Saito, Masami Y. Hirai, and Keiko Yonekura-Sakakibara. “Decoding genes with coexpression networks and metabolomics–‘majority report by precogs’.” Trends in plant science 13.1 (2008): 36-43. [42] Joshua M. Stuart, et al. “A gene-coexpression network for global discovery of conserved genetic modules.” Science 302.5643 (2003): 249-255. [43] Claire M.M. Gachon, et al. “Transcriptional co-regulation of secondary metabolism enzymes in Arabidopsis: functional and evolutionary implications.” Plant molecular biology 58.2 (2005): 229-245. [44] Keiko Yonekura-Sakakibara, et al. “Identification of a flavonol 7-Orhamnosyltransferase gene determining flavonoid pattern in Arabidopsis by transcriptome coexpression analysis and reverse genetics.” Journal of Biological Chemistry 282.20 (2007): 14932-14941. [45] Koh Aoki, Yoshiyuki Ogata, and Daisuke Shibata. “Approaches for extracting practical information from gene co-expression networks in plant biology.” Plant and cell physiology 48.3 (2007): 381-390. [46] Harold W. Cook. “Fatty acid desaturation and chain elongation in eucaryotes.” New Comprehensive Biochemistry 20 (1991): 141-169. [47] Kim D Pruitt, et al. “RefSeq: an update on mammalian reference se- 61 quences.” Nucleic acids research 42.D1 (2014): D756-D763. [48] Claire R. Shen, et al. “Driving forces enable high-titer anaerobic 1-butanol synthesis in Escherichia Coli.” Applied and environmental microbiology 77.9 (2011): 2905-2915. 62 発表リスト ポスター発表 [1] 奥田裕樹, 小野直亮, 有田正規, 太田大策, 佐藤哲大, 杉浦忠男, Md.Altaf-UlAmin, 金谷重彦, “非モデル生物における遺伝子発現解析のための de novo トラン スクリプトームアセンブリの性能評価”, NGS 現場の会 第三回研究会, 2013 年 9 月 4,5 日 [2] 奥田裕樹, 小野直亮, 有田正規, 太田大策, 佐藤哲大, 杉浦忠男, Md.Altaf-UlAmin, 金谷重彦,“ミドリムシにおける次世代シーケンサを用いた遺伝子発現差異解 析”, 第 6 回生命情報科学若手の会 研究会, 2014 年 10 月 29,30 日 63 付表 表 4.1: 藻類データベースに含まれる生物種名と配列数 Species Name # of sequences Acutodesmus obliquus 97 Aureococcus anophagefferens 105 Aureoumbra lagunensis 110 Bryopsis hypnoides 69 Chaetosphaeridium globosum 144 Chara vulgaris 151 Chattonella marina 41 Chlamydomonas moewusii 14 Chlamydomonas reinhardtii 14,488 Chlorella variabilis 80 Chlorella vulgaris 174 Chlorokybus atmophyticus 172 Chrysodidymus synuroideus 37 Coccomyxa subellipsoidea C-169 111 Cryptomonas paramecium 548 Desmarestia viridis 39 Dictyota dichotoma 38 Dunaliella salina 91 Durinskia baltica 129 Ectocarpus siliculosus 148 Emiliania huxleyi 140 Fistulifera sp. JPCC DA0580 135 Floydiella terrestris 74 Fucus vesiculosus 177 Gonium pectorale 76 64 表 4.1: 藻類データベースに含まれる生物種名と配列数 Species Name # of sequences Guillardia theta 632 Helicosporidium sp. ex Simulium jonesi 63 Hemiselmis andersenii 515 Heterosigma akashiwo 197 Kryptoperidinium foliaceum 139 Laminaria digitata 39 Leptosira terrestris 88 Mesostigma viride 146 Micromonas pusilla CCMP1545 10,269 Micromonas sp. RCC299 10,140 Monomastix sp. OKE-1 82 Nannochloropsis gaditana 152 Nephroselmis olivacea 195 Nitella hyalina 48 Ochromonas danica 44 Oedogonium cardiacum 99 Oltmannsiellopsis viridis 129 Ostreococcus lucimarinus CCE9901 Ostreococcus tauri 7,603 398 Parachlorella kessleri 84 Pavlova lutheri 111 Pedinomonas minor 93 Phaeocystis antarctica 100 Phaeodactylum tricornutum 165 Phaeodactylum tricornutum CCAP 1055/1 10,408 Polytomella capuana 7 Polytomella parva 7 65 表 4.1: 藻類データベースに含まれる生物種名と配列数 Species Name # of sequences Polytomella sp. SAG 63-10 7 Prototheca wickerhamii 36 Pseudendoclonium akinetum 177 Pycnococcus provasolii 86 Pylaiella littoralis 52 Pyramimonas parkeae 94 Rhodomonas salina 190 Saccharina angustata 38 Saccharina coriacea 38 Saccharina diabolica 38 Saccharina japonica 177 Saccharina japonica x latissima 38 Saccharina longipedalis 38 Saccharina ochotensis 38 Saccharina religiosa 38 Schizomeris leibleinii 77 Staurastrum punctulatum 103 Stigeoclonium helveticum 79 Synedra acus 165 Thalassiosira oceanica CCMP1005 142 Thalassiosira pseudonana 176 Thalassiosira pseudonana CCMP1335 Vaucheria litorea 11,673 139 Volvox carteri f. nagariensis 14,436 Zygnema circumcarinatum 103 66