8月びわこ その2

　今日は今月二回目の琵琶湖。ここ最近台風や豪雨もなく、ずっと猛暑が続いていたので、あらゆるものが「真夏の琵琶湖の最強版」になっていた。

　まず船着き場についてすぐに目についたのがアオコ。水の動きが少ないと発生しやすいと言われているけど、ずっと雨が降ってないせいで、琵琶湖の水位も-36 cmまで下がっていて、南郷洗堰の放水量もずっと30トン以下。この夏、琵琶湖の水はほとんど動いていない。昨日見てきたけど、瀬田川までもアオコで緑の水になってしまっている状況だ。

　その一方で、北湖では水がますます透明になっている。少雨で川からの泥の粒子の流れ込みがないうえ、猛暑で温度躍層の発達が進んで、表層の栄養不足による植物プランクトンの減少がすすんでいるからだ。水深10 mまでが水温28℃以上、そこから水温が一気に下がって、水深20 mではもう13℃という状況。透明度は8.5 mで、採水器をおろしても外洋のように水が透き通っている。

　そして水が透明になると現れるのが最小の植物プランクトンであるピコシアノバクテリアだ。5 mで採水した水を早速蛍光顕微鏡で見てみたのだけど、うじゃうじゃいた。

ちなみに同じ視野でみた細菌がこれ↓

　細菌が4.6×10^6 cells/mlに対して、ピコシアノバクテリアが3.0×10^5 cells/ml もいて、全細菌の6.5%を占めている計算。

　ピコシアノバクテリアが元気に増えているということは、それだけそれに感染するウイルスもたくさん発生しているはずだ。ということでウイルスを見てみると、これまた今までに見たことが無いくらいたくさんいた(弱く光っているやつも結構います)↓

　現存量は1.2×10^8 particles/mlで、気のせいかもしれないけど、粒子の光り方もいつもより強くて、状態が良いものが多いような気がした。真夏の透明な水、面白い。透明ということは、大きな生き物や川から流れてきた泥の粒子が少ないということでしかない。小さな生物やウイルスの世界は、透明な水の中でむしろ活発に展開されている。今日採ったサンプルのシーケンスデータの解析が楽しみだ。

　そして最後は積乱雲の発生で急に波風強くなり、慌てて撤退。何から何まで真夏でした。

　次の日曜日からは海外出張。カナダで行われるISMEで発表と情報収集してきます。すでに要旨が発表されているけど、僕が研究しているマニアックな細菌を研究している人がいることが判明して、しかも僕が研究している環境とは違う環境から見つかったという内容らしいので、どんな発見が発表されるのか、すごく楽しみだ。

水圏微生物研究フォーラム

「水圏微生物研究フォーラム」なる集会に参加するため、関東出張してました。

　水域の環境微生物の研究者のための集会で、昨夏に参加した国際学会のSAME (Symposium on Aquatic Microbial Ecology)の国内版という感じ。招待講演のメンバーを見た瞬間に「こんなん絶対おもろいわ」と思って即座に申し込みをしたのだけど、期待通りで、ポスター含め、一つたりとも聞き逃すわけにはいかない濃密な内容だった。あまりに面白いポスターばっかりだったので、時間外もずっとポスターを見て回っていたのだけど、それでもちゃんと聞けずに終わってしまったポスターがたくさんあって、もっとポスターの時間を伸ばしてほしいと思ったほどだ。

　自分のポスターを聞きに来てくれる人も近い分野の人が多いので、密な議論ができて、「この話で盛り上がれる日本人はみんなここにいるのでは？」というようなレベルのマニアックな会話もできたのはとても楽しかった。

　またよく言われることだけど、少人数の学会って、近い分野の研究者と密に議論できるというだけでなく、大規模学会ではなかなか話すタイミングが得られない、大御所の先生方と話すチャンスも結構あるというメリットもある。僕自身も今回、是非お話ししたいと思っていた研究者とお話しすることができて、とても満足だった。

　総じて、これまでに参加した国内の研究集会では1番良かったと言ってもいいくらいの満足度だったと思う。規模的にも、時間的にも、ちょうど良い内容だった。でかい学会で異分野の研究を見るのも楽しいし、業界のトレンドを追うのには良いのだけど、やっぱりこういう「狭くて深い」研究会のほうが、議論の質も、得られる人脈の質も上だと感じた。今後もこういうチャンスがあったら是非参加したい。

8月びわこ その１

今月は2回琵琶湖に行く予定なので、早速1回目を済ませてきた。

穏やかな夏の朝・・・かと思ったけど、やっぱり日が昇ると暑い。9時の時点で気温は29.7℃、水温は28.1℃。透明度は7.9mまで伸びて、「溶存有機炭素（DOC）豊富、栄養塩（N・P）欠乏」の典型的な真夏の表層になってきた。

沖ではカモメが表層の魚を追いかけていた。今年の琵琶湖はアユが多いらしいので、おそらくアユの群れだと思う。その影響か、今年は水中のミジンコの数が去年に比べて遙かに少なく、去年は透明なボトルに水を汲むと泳いでいるのがたくさん見えたのに、今年はそんなに見つからない。南湖の水草も例年のように水面を覆いつくすほどではないし、烏丸半島のハスが突然消えたのもニュースになった。一方で外来植物のオオバナミズキンバイの黄色い花が去年にも増してそこら中の岸に生えているし、去年は見られなかった家の近所の瀬田川でも今年は大繁殖している。

1年単位でも、琵琶湖はどんどん変わっている。変化のスピードは人間社会とそんなに大差ないのではないかとも思ってしまう。こんな勢いで数百万年も続いている琵琶湖とその変化に順応して生息し続けている生物って本当にたくましい。

初稿という既成事実

　月例の調査やサンプル処理は進めつつも、最近はデスクワーク中心の毎日。次の論文に載せるデータを作りながら、そのデータを使って来月の国際学会での発表ポスターを作成中。

　前回の論文はかなりスピード重視で書いたけど、今回のはあえて焦らずに分析を深める方向で進めている。とにかく早くアウトプットが欲しかった前回と違って、今回はスピードにそこまで貪欲ではないということや、前回よりもデータ量もストーリーもでかいので分析に時間がかかるということもある。だけどそれ以上に、前回の論文執筆・査読対応を通じて

丁寧な分析というのは思ったよりも評価されるものであり、初稿の段階で分析を深めておくことがむしろ、労力に対して質の高いアウトプットを出すことにつながるのではないか

という思いを強くしたことが理由としては大きい。

　前回の論文は、会社員時代の仕事の仕方を引きずっていたこともあって、

判断に迷うところはとりあえず仮決めで進めて、早めに審判を仰ぎ、軌道修正が必要かどうかを早く知るほうが効率がいい

という考えで、あまり深い方向に突っ走らず、できるだけ最小公倍数的なモノを書いて査読に回すことを目指してすすめていた。その結果、根本の考え方のエラーを指摘され、分析まるごとやり直しになって、二度手間を回避することができたパートもあった。だけど同時に面白かったのが、自分は「最小公倍数で仮決めだ」と思っていた分析でも、評価されて結局最後まで生き残ったパートも結構たくさんあった、ということだ。必要であればもっと深い分析に突っ走ることもできたけど、そこまで行く前に出しておいて意外に大丈夫なんだな、ということを思った。もちろん、掲載先のジャーナルのレベルにも大きく依るのだろうけど、結果として「とにかく早く出したい」という僕の欲望に応える結末になった。

　これは決して「手抜きで不完全な論文を出してしまった」ということではない。言いたいのは「もっと深い方向に突っ走ることもできたけど、敢えて深入りせずに置いておいたら、それで意外とうまくいった」ということだ。仮に深い分析をして出していたところで、「不要なこだわり」と判定されて、ボツになっていた可能性だって大いにあるだろう。だけど一通り論文が出てみて「最初からできる限り深い分析をしておいても良かったな」というのもちょっと思った。

会社にいた時も思っていたけど、

「仮決めで仕事を進めること」で二度手間を回避するチャンスも増えるけど、誰かが修正してくれるだろうと思ったら誰も指摘してくれずにそれがそのまま既成事実になってしまうリスクも増える

ということが、研究論文にも当てはまると感じた。いくら業界のプロが査読するとはいえ、「査読⇔修正」のやり取りはせいぜい数回までだし、時間が限られている中で、初稿がベースになって議論が進んでいくのはまあ当たり前のことだ。

　だから今回は、最初から分析できるところまで分析しつくした、最高クオリティの初稿を書き上げたいと考えている。判断に迷うところも、仮決めで進めるのではなく、全ての可能性を盛り込んだ最大公約数的な分析にしていきたい。その結果、出したものの多くが不要になるかもしれないけれど、深く分析したこと自体は査読者はちゃんと評価してくれて、無駄になることはなさそうだという自信が出てきたことと、初稿で最高クオリティのものを送っておかないと、後からその既成事実を上回るものに書き換えることは難しいということが分かったからだ。

　なので、当初の予定からかなり後ろ倒しになっているけど、緻密な議論・分析をすべく、じっくりと時間をかけて次の作品の準備を進めている。「通ればいい」という欲望ではなく「驚きとワクワクを与えたい」という欲望をモチベーションにして、恥ずかしくないものを世の中に出したいと思う。

7月びわこ

今日は琵琶湖調査の日でした。ちょっと波風があるけど、良い天気。

気温28度でまだマシな暑さ。水温は26度、透明度は6mに増えて、水中も真夏モードに入りつつある。

いつもは深層ばかり採水しているのだけど、今回は別の目的で久しぶりに表層から採水。やっぱり深層と全然違って粒子が多いので、処理の大変さが全然違う。いつもの半分の濾過量でフィルターが水漏れを起こしてしまった。

さっそく細菌数を数えてみたのだけど、7.1×10^6 cells/mlで、深層の7倍近い現存量。これは琵琶湖にしてはかなり多い数だけど、盛夏に向かっていく今の時期には過去にも報告されているレンジでもある。梅雨が終わって、表層に残された最後の栄養の取り合いで細菌達も数を増やすのに必死なのかもしれない。ウイルスは4.8×10^7 particle/mlで、これも深層の倍以上の値。

変化が激しい表層の調査は色々緊張する。深層の安定性や再現性があってロバストなデータが出る感じは、やっぱりいいよなぁと思った。

琵琶湖に生息する微生物の図鑑

自分の研究の論文が出版されたので、簡単に紹介。

Vertical partitioning of freshwater bacterioplankton community in a deep mesotrophic lake with a fully oxygenated hypolimnion (Lake Biwa, Japan) - Okazaki - 2016 - Environmental Microbiology Reports - Wiley Online Library

一言でいえば

琵琶湖の水中で季節・深さごとに出現する細菌を網羅した図鑑を作成

3行で言えば

・琵琶湖沖で15か月・3水深にわたって、水中の細菌群集を網羅的に調査

・特に、研究例が少ない、深層に生息する細菌に注目

・深層にのみ生息する、10以上の「変わった系統」の細菌を特定

という内容。

研究する意味

「琵琶湖の細菌の研究なんてもうやりつくされているんじゃないの？」と思われそうだけど、実はまだほとんど何もわかっていない。単に研究者がいなかった、というのも大きな理由だけど、やはり一番大きな障壁は、 

ほとんどの細菌は培養ができないので、研究もできない

ということだ。

琵琶湖の水中には1ml（ティースプーン1杯）あたり100万を超える細菌が生息している。顕微鏡でみれば、それこそ星の数ほどの細菌が見られる。

　だけど、環境中にはこんなにたくさんいるのに、なぜか研究室に持ち帰って人工的な環境で飼おうとすると、ほとんどの細菌は増えてくれない（正確には、少数の増えやすい奴ばかり増えてしまう）。その理由はまだよく分かっていない。人工的な環境では再現できないような、環境中に存在する絶妙な何かが欠けているのだろう。で、培養ができなければ、個々の細菌系統の性質を詳しく調べることもできない。だから、難培養の細菌系統を培養する、というのは大仕事で、日夜多くの研究者たちがリソースを割いているし、重要な系統を単離することに成功すれば、それだけでNatureに掲載されるレベルのインパクトを世界に与えることができる。

　琵琶湖の細菌も例外ではなく、細菌全体の存在数やその季節変動を報告した研究は昔から例があったものの、培養できないことが障壁となって、「どのような系統が生息しているか」ということは解明されずにいた。

　そこに訪れた大きな転機が「次世代シーケンサー（NGS）」の登場だ。簡単に言えば、環境中に存在するDNAを集めてきて、片っ端から塩基配列を読むことで、どの遺伝子を持った細菌がどれくらいいるかを、培養することなしに推定することができる技術だ。細菌の系統は特定の遺伝子（16S rRNA遺伝子）の塩基配列を元に分類でき、データベースに情報が蓄積されている。だから、この遺伝子の塩基配列の情報さえ手に入れば、それが培養できない系統であっても、どのような分類に属し、過去にどのような環境で見つかってきたものなのかは分かる。

 このNGSの登場は本当に革命的で、水中に限らず、あらゆる環境に生息する細菌の生態が急速に明らかになりつつある。微生物学者にとっては「ボーナスステージ」と言っていいような状況で、様々な動物の腸内、家のドアや電車のつり革、歴史的絵画の表面に至るまで、調べれば調べるだけ新しい発見が出てくる事態で、論文もバンバン出版されている状況だ。

 その一方で、「NGSを使ってまだ調べられていない環境の細菌の組成を調べました、おしまい」というだけの、単なる記載的な研究が増えているという批判があるのも事実だ。

新しいところを調べれば、新しいものが出てくるのは当たり前。それが何を意味しているのか、どう面白いのかが肝要だ

という流れになりつつあるというのが、NGS技術が一通り普及しきった今の業界の状況だ。

何をしたのか？

　このNGSを琵琶湖で初めて使い、どのような細菌が生息するのかを網羅したのがこの研究だ。だけど、単に表層の水を採って細菌群集を調べるだけの研究であれば、すでに世界中の湖沼で行われている研究の二番煎じになってしまう。

そんな中で、本研究は、

①これまでに研究が行われてこなかった、深層に焦点を当てたこと

②15か月にわたる継続的な調査を行い、季節変動を明らかにしたこと

という二点を大きなウリにしている。

特に①の「大型淡水湖の深層にはどのような細菌が生息しているのか」ということに関しては、世界的にもほとんど報告が無い。手掛かりになるような断片的な情報は存在するものの、NGSで季節変化を含めて網羅した研究は、本研究が初めてになる。そのような意味で、本研究は

湖の深層に生息する細菌の組成を、世界で初めて網羅的にみる 

のが目的と言える。

何が分かったのか？

　過去に琵琶湖沖に毎月通って採り溜めていた3水深×15か月のDNAサンプルを使って、NGSで16S rRNA遺伝子のシーケンスを行い、とれてきた塩基配列をもとに各サンプルの細菌群集を分析した。すると、表層と深層では生息する細菌系統が大きく異なることが分かった。しかも、これまで表層で良く知られている細菌とは、「門」のレベル（一番大きな生物の分類の区分）で異なる系統が深層にたくさん生息していることが分かった。

　これらの系統の塩基配列をDNAのデータベースで調べてみると、実はすでに、世界中の深い湖からの報告が散発的に存在することが分かった。このことから、これらの系統は、琵琶湖にだけ生息するのではなく、世界中の湖の深水層に普遍的に生息する系統であると考えられる。

　本研究では、これらの過去に報告されてきた知見をまとめながら、琵琶湖で季節・深度ごとに採られたデータを加えて、湖沼の深水層に生息する細菌系統を初めて網羅的に整理することができた。

何がすごいのか？今後の発展は？

　湖の深層に生息する細菌は、湖全体の物質循環や生態系の中で重要な役割を果たしている。本研究が「深層には表層とは全く違う系統が生息している」ことを明らかにしたことで、今後これらの深層独自の細菌系統の生態の解明を進めていかなければならないことが分かった。

　ただ、「じゃあ彼らは生態系内でどのような役割を果たしているのか？」とことについては、現時点では分からない。今回見つかった「門」の多くは、多数の難培養性系統で構成される、ほとんど正体が分かっていない分類群だからだ。残念ながら本研究は「変なのがたくさんいるからもっと研究を進めましょう」というところで終わってしまっている。

　なので、今後は今回明らかにした個々の細菌系統に焦点を当て、どのような遺伝子や代謝活性を持っているか調べて、彼らが生態系で何をしているのか、なぜそこにたくさんいるのか、ということを明らかにしていく必要がある。

　また、今回は琵琶湖のみでの調査だったけど、データベースの情報などから、他の湖にも似たような系統が存在している可能性が高い。今回琵琶湖で見つかった発見が、どの程度一般化できる内容なのか、今後調査地の数を増やして検討していく必要がある。

　前者については、データをとり始めており、後者についてはすでに結果が出揃って、論文の執筆に着手しているところだ。

　また、さらに長い目で見たとき、僕はこの湖の深層の細菌の研究を、単なる湖の生態学の研究で終わらせずに、「生物学の新発見の宝庫」に育て上げたいと思っている。断言できる根拠はないけれど、環境の特殊性や、調査の容易さ、出現する細菌系統の異質さなどを考えると、この環境を深堀りしまくることで、生態学や微生物学、もしかすると生命科学全般に影響を及ぼすような、面白い現象を発見できるのではないか、という気がしている。本研究は、人生を賭けるその長い道のりの、最初の出発点となる予定の研究だ。

ここに至るまで

　しかし、長かった。会社を辞めて研究に戻ってきてから、1年以上アウトプットが無かったことがとてもストレスだった。ライバルが似たような研究をしているのも知っていたし、何よりも「色々やっているのにまだ何も成果になっていない」ということが悔しくて、この1年少し、早く論文を出したい一心だった。

　本研究は、去年の4月、研究に戻ってきた直後に思いついて、修士課程の時に採り溜めていたサンプルを使ってすぐに着手した研究だ。「とりあえず論文1本出そう」くらいの気持ちで、そこからほぼ滞りなく、最速で進めてきたつもりだけど、

・4-5月　予備実験

・6-7月　シーケンス外注

・8-9月　結果分析、一部サンプルの失敗判明　　

・10-11月シーケンス再外注

・12月　　分析

・1月　　論文執筆⇒英文校閲

・2月　　論文投稿⇒瞬殺リジェクト⇒すぐに再投稿

・3月　　レビューコメント付きリジェクト（姉妹誌への投稿はOK）

・4月　　リバイス⇒姉妹誌に投稿

・6月　　マイナーリビジョン⇒すぐに再投稿⇒アクセプト

・7月　　オンライン公開

　という流れで、結局1年以上もかかってしまった。特に、初稿を書き終えた1月から、アクセプトまでに半年も待たされたことは、とてつもなくストレスだった（それでもかなり早いほうだと言われるけど）。また、2回目のリジェクトのあと、査読コメントに従って大幅に書き直しをすることになった4月は本当に苦しかった。一度書いたものを全て分析しなおして、書き直す大変さは並大抵ではない。初稿ができた段階で、8割くらい終わっているかと思ったけど、時間的にも労力的にも、初稿ができた段階でようやく半分くらいだなぁ、と思った。

　会社員時代と同じ仕事スピードで大学院生をやれば、サクサク早く効率的にアウトプットが出せるだろうと思っていた頃もあった。だけどやっぱ研究は、厳密さが全然違うし、一人でやらなければならないから、必要な労力も時間も全然違う。種を撒いてから収穫まで、1年以上は余裕でかかる。だから、去年の秋に種を撒いて、今論文を書いているネタが公開されるのは来年になるだろうし、今種を撒いているやつが出荷されるのは再来年の春とかになるかもしれない。本当に気が遠くなるし、のろまだと思う。でも、しょうがない。ライバルもみんな同じだ。効率的に収穫できるように、できるだけ計画的に種まきをすすめるのみだ。

　しかし、やっぱり、自分の仕事が自分の名前で論文になるのはとても嬉しい瞬間だ。一本目の論文が出た時も「これでいつ死んでも世の中に名前が残せる」みたいなことを思ったけど、またこの感動を味わえて嬉しい。研究に戻ってきてよかったです。

RによるNMDSを用いた微生物群集構造解析

今日も解析で苦労したので備忘メモ。
※今後、自分の理解に合わせて勝手に加筆・修正するかもしれません。

以下のものは「とりあえず動かすところまで」を目標に書いたものです。もし間違いがあった場合はご指摘くださると有難いです。

なお、今回勉強するにあたっては、以下の記事を大いに参考にしました。英語がOKであれば、こっちを読んだ方が正確で手っ取り早いかと思います。
NMDS Tutorial in R – sample(ECOLOGY)
NON-METRIC MULTIDIMENSIONAL SCALING (MDS)
Ordination analysis [Analysis of community ecology data in R]

目的

　Rを使って、16Sアンプリコンシーケンスの分析から得られたOTUテーブル（site = 50 / OTU = 800 程度）から、NMDS(non metric multidimensional scaling)を用いてサイト間の群集の類似性（β多様性）を比較し、さらに各サイトの環境要因のメタデータ（温度・栄養塩・酸素濃度・・・）の情報を載せて、どの要因がサイト間の群集の差を強く説明しているのかを一覧できる散布図を作る。

こんな図を作りたい（論文でよく見るやつ）↓

なぜNMDS？

　複雑なものを2次元のプロットに落として視覚化するタイプの多変量解析にはPCAやCAなどの手法もある（参考※pdfです）。NMDSは微生物生態学の論文では良く見るけど、他の分析と比べて一体何がいいのだろうか？
　一番の利点は、rankに基づいて距離を計算しその距離関係に基づいてプロットを描くというノンパラメトリックな手法であるため、より幅広い事例に応用できてロバストな結果が出やすい、ということだ。さらに目的に応じて距離の計算方法を自由に変えられる、ということも利点とされている。（後者に関しては、以下で紹介する手法のデフォルト設定値（Bray-Curtis指数による計算）がすでに生態学ではほとんどの場合最適とされる計算方法になっているため、今回は触れない（とくに設定を触る必要は無い）。）
　逆に欠点としては、NMDSが反復的に計算して最適解を見つけてくる手法であるため、他の手法のように1点の最適解に落ち着かず、計算の度に結果が微妙に変わることや、マシンパワーを必要とすることが挙げられる。これに関しては、以下で示すように、ランダム試行の出発点となるseed値と、計算の反復回数を固定することで、再現性を担保した結果が得られるようにすることはできる。また各軸がどれだけデータのばらつきに寄与しているかが示されるPCAなどと違って、NMDSでは軸自体には何の意味も無いため、点が近いか遠いかの位置関係だけによって結果が解釈されるのも注意点だ。

準備するもの

Rにパッケージ"vegan"を入れておく。
以下2つのデータを行名・列名を付けた状態でデータフレームにしてRに読み込ませる。

• ①OTU_table (各行がsite,各列がOTU)
• ②Environmental_parameters (各行がsite,各列が環境要因)

Rarefactionを行う場合

　通常アンプリコンシーケンスを行うときは、各siteのDNA量が均等になるようにサンプルを準備するため、本来はどのsiteも同じ数だけDNAリードが取れてくるべきなのだけど、このゆらぎが結構大きいようで、往々にしてサンプル間のリード数が最大3倍くらい違うシーケンス結果が出てくる。リード数が違うまま異なるsiteのサンプル同士を比較すると、当然バイアスがかかる。なので、そのばらつきを平準化するために、各siteからリードをランダムにリサンプリングするのがRarefactionだ。この手法、微生物生態学では広く行われてきた作業だけど、得られたリードの一部しか使わない（往々にして半分以上捨てることになる）、言い換えると（ランダムに選ぶとはいえ）「とったデータを無かったことにする行為」になるため、不正確な分析につながっているという問題点も指摘されていたりするので、（特にレアなものを対象にしたい場合など）データの扱いには慎重になる必要がある。
　なお、古くから一般的に行われているRarefactionは、「最もDNAリード数が少ないサンプルに合わせて他のサンプルのリード数を削る」というAbundanceベースのRarefactionなのだけど、最近になって「母集団の多様性推計値をどれくらいカバーできているかで各サンプルのリード数を合わせる」というCoverageベースのRarefactionも提唱されていて、こちらのほうがサンプル間の多様性の違いを正確に比較できるという報告もある（文献）。まずは両方やってみる、ということで、以下、それぞれの方法を説明。

AbundanceベースのRarefaction

　veganに入っているrrarefyという関数でランダムリサンプリングができるので、リード数が最も少ないサンプルに合わせてサンプリング数を指定してやれば良い。

set.seed(123) #再現性をとるためにランダム変数を固定（数字は何でもいい）
OTU_aburared<-rrarefy(OTU_table,min(rowSums(OTU_table))) #最小サンプルに合わせてリサンプリング
write.csv(OTU_aburared, "abundance_rared.csv") #CSVに出力したい場合

CoverageベースのRarefaction

※以下に関してはメタバーコーディング・環境DNAバーコーディング解析の技法の「４．メタゲノムデータを用いた群集統計解析法」の掲載資料を大いに参考にさせていただきました。

　Coverageベースでやる場合、最後は同じくrrarefyでリサンプリングをしてくるのだけど、その前にRarefactionカーブの傾き（＝カバー率）が一定の値（＝一番カバー率が低いサンプルの傾きと同じ傾き）以下になるDNAリード数をサンプルごとに求めなくてはならない。veganのrareslopeという関数を用いて、Rarefactionカーブの傾きを求められるので、各サンプル、1リード増やすごとの傾きを、forを使ってとってくるところがポイント。

#各サンプル、１リード増えるごとの傾きをリスト形式で一括出力（要計算時間）
rareslopelist<-list()
for(i in 1:nrow(OTU_table)){
  rareslopelist[[i]]<-rareslope(OTU_table[i,],1:(sum(OTU_table[i,])-1))
}#最後の1リードを入れずに「-1」としているのは、最後の1リードについては傾きが計算できない（0になる）ため

#最後のリードの傾きが最も大きい（＝カバー率が最も低い）サンプルを探す
getmincov<-c()
  for(i in 1:nrow(OTU_table)){
    getmincov[i]<-rareslopelist[[i]][length(rareslopelist[[i]])]
}#最後の1リードの傾きだけをベクトルでとってくる

(1-max(getmincov))*100　#カバレッジの％を確認

#指定したカバレッジに到達した（＝傾きが指定値を下回る）瞬間のリード数をサンプルごとに採ってくる
cvrfun<-function(x){min(which(x<=max(getmincov)))+1} #関数を設定。上記で1を引いた分を足し戻す
cvrrare<-unlist(lapply(rareslopelist,cvrfun))　#lapply+unlistでベクトル形式にして一括で値を取得

set.seed(123) #再現性をとるためにランダム変数を固定（数字は何でもいい）
OTU_covrared<-rrarefy(OTU_table,cvrrare) #得られたリード数に沿って各サンプルからリサンプリング

write.csv(OTU_covrared,"covarage_rared.csv")  #CSVに出力したい場合

（補足）α多様性の確認

diversity関数で各siteの群集そのものの多様性（α多様性）の計算も可能。変数を指定することで、複数の多様性指数(ShannonとかSimpsonとか)を計算してくれる。

diversity(OTU_aburared,index="shannon")

NMDSをかける

分析にかけるデータフレームができたら、"vegan"のmetaMDS関数を使って、NMDSプロットに必要な「距離の計算」を行う。

set.seed(123) #NMDSの再現性を担保するため、ランダム変数を固定（数字はなんでもいい）
nmds<-metaMDS(abundance_rared,k=2,trymax=20) 
#"abundance_rared"の部分はphyseq-classのデータでもOK
#kは次元数。基本は2次元にしたいので2を使う
#trymaxは計算の最大試行回数。増やすとより良い解が見つかりやすくなるが、
#計算時間もかかるため、検討段階では少なめに、本番では大きな値を使ってやる
#この回数に達する前にプログラムが最適と判断して計算を止めることもある

ここで、NMDSの出来を評価するのが"stress"という値だ。これが小さければ小さいほど、データを上手く2次元に落とし込めたということを意味しているらしい（原理は完全には理解していないので不正確な表現かもしれません）。この値は、

nmds$stress

で見ることができる。で、どの程度であればOKか？という話なのだけど、Web上の情報で良く出てくるのは「>0.3だとダメダメ、>0.2だとまあまあ、>0.1だと良い、>0.05だと最高！」という感じの指標なのだけど、サンプルサイズや、site間の環境の違いの大きさなどにも影響するため、一概には言うべきではない、という主張もあったりして、統一見解はなさそう。ので、できるだけ低いstress値が出た結果を、その値と一緒に載せるのが外に出すときの一般的なやり方っぽい。
　なお、set.seedの値をいじれば分かるけど、最初の出発点によって結構結果が変わってくる。出発点によっては局所最適みたいなところに落ち着いてしまって、計算回数をいくら増やしてもstressが減らないということも起こる。ので、本番の分析をやるときは、

stress.value<-NA
for(i in 1:100){
　set.seed(i)
　nmds<-metaMDS(abundance_rared,k=2,trymax=50)
　stress.value[i]<-nmds$stress
}
grep(min(stress.value),stress.value)

とかして、低いStress値を叩きだすseed値を探してくるのが良さそう。
　ここまでくれば、あとは作図関数に投げるだけだ。veganに入っているordiplot、orditorpを使えば、metaMDSで生成したデータ（上記でいう"nmds"）を、縮尺など調整して綺麗にプロットしてくれる。

win.graph(30,30) #正方形の作画領域を生成
ordiplot(nmds,type="n") #枠を書く
orditorp(nmds,display="sites",air=0.1,cex=1)　#点を書く
#displayをsitesとすればsiteの、speciesとすればOTUのプロットを描いてくれる
#airは重なって表示されるラベルをずらすための変数
#変数labelsで表示ラベルを、selectで表示/非表示プロットを指定可能
#普通のplotの変数がそのまま使え、好きな色や形の点を打てる

環境変数をNMDS上に載せる

　ここまでで各site間の組成がどの程度違うのか（β多様性）を、視覚化することができた。近いものほど似ており、離れているものほど違う。この段階で誰の目にもはっきり分かるくらい傾向が出ていれば、これ以上の説明はいらないだろう。だけど、一目で結果が分かりにくい時、各siteの環境条件のデータを付け加えることで、「温度が強く効いている」とか「酸素濃度はそれほど影響がなさそう」とかが見えてくることがある。そこで冒頭に示したもう一つのデータセット「②Environmental_parameters」をveganのenvfit関数を使ってNMDS上に反映させる。

envNMDS<-envfit(nmds,Environmental_parameters)
#データにNAがある場合は変数"na.rm=T"を加える

これで先ほど作った"nmds"上に"Environmental_parameters"の値を載せた場合の各環境変数の説明力が計算された。このデータを見れば、各環境変数が、どの程度NMDS上の各siteの距離関係のばらつきを上手く説明できているかを示す表が得られる。

envNMDS
#             NMDS1    NMDS2     r2 Pr(>r)    
#env1       0.77127 -0.63650 0.4564  0.001 ***
#env2      -0.67797  0.73509 0.6592  0.001 ***
#env3      -0.07703 -0.99703 0.1305  0.175    
#env4      -0.15870  0.98733 0.5373  0.001 ***
#env5       0.56661  0.82399 0.3727  0.004 ** 

当てはまりの強さは、P値で示されており、上の例ではenv1,env2,env4がうまくsite間のばらつきを説明できていて、env3が説明できていない環境要因だということが分かる。
これを先ほどのNMDS上にプロットする。普通のplot関数で"envNMDS"を読んで上記のNMDS上に上書きするだけで、矢印まで描いてくれるようになっている。

win.graph(30,30) 
ordiplot(nmds,type="n") 
orditorp(nmds,display="sites",air=0.1,cex=1)　
#ここまでは上記と同じ
plot(envNMDS,p.max=0.01) #p.maxでその値以下のP値を叩きだした環境要因のみプロットする

これで冒頭の図と同じものができる。なお、各環境要因の矢印だけど、

向きと相対的な長さだけが意味を持ち、矢印の始点と終点、絶対的な長さには意味が無い

ということに注意。まず矢印の向きだけど、「この向きがsite間のその環境要因の傾きが最も急な方向である」ということを示している。そして矢印の長さは「その傾きの急さ」を示しており、上記のr2の値の平方根に比例してスケーリングされるようになっている。そしてこの「長さ」は各矢印間の相対値で見るべきで、絶対的な長さは、プロット上で最も見栄えが良いように自動的に決定されているに過ぎない。つまり、矢印はNMDS上のプロットの縮尺とは独立した縮尺で描かれていて、その位置関係を比較することはできないということだ。
　この説明では少しわかりにくいかもしれないけど、同じくveganについているordisurfという関数を用いることで、この傾きの「傾斜」と「向き」を等高線的に表示することもできて、これを見れば矢印がどういう意味なのかももう少し分かり易いのではないかと思う。

ordisurf(nmds,Environmental_parameters$env1,add=T)
#"add=T"によって、先述ののNMDSプロットに上書きできる。
#インプットデータは元のデータフレームから１列のみ抜き出したものにする

この関数の原理については、僕はきちんと理解できていないのだけど、以下のページが大変詳しい。
What is ordisurf() doing...?
　なお、これらの情報を調べていると「ノンパラメトリックな手法であるNMDS上にenvfitで出したリニアな回帰の結果を載せていいのか？」という議論を所々で見かけるのだけど、veganの開発元でも本件はFAQとして取り上げられていて、「距離の計算方法はノンパラメトリックであるが、NMDSのプロット自体は距離スコアに応じた線形のデータなので、envfitやordisurfで作ったデータを載せることは問題ない」という見解のよう。一方で、この方法で論文のレフェリーに突っ込まれたという話も出てきたりして、自分の直感を超えた統計手法を使うときは、目的をはっきりさせて原理としっかり理解してからでないと怖いなぁということで、まだまだ勉強途上ですが、ここまでの理解をまとめてみた次第です。