両生類の個体群変動の大きさに影響を与える要因


草野　保

　はじめに

　近年、地球規模での両生類の減少が危惧され、日本でも、特に大都市近郊でのこの数十年の両生類個体群の絶滅・衰退は際立っている。実際、環境省のレッドリストでは、有尾類12種、無尾類10種が絶滅のおそれがある種（絶滅危惧II類以上）としてリストアップされている（環境省, 2014）。現在国内には、有尾類35種と無尾類43種の在来種が生息するので（日本爬虫両棲類学会, 2018）、絶滅危惧種の割合はそれぞれ34と23％となり、特にサンショウウオ類は著しく高率となっている。両生類の絶滅や減少の原因には、人為的な環境破壊に伴う生息環境の喪失とともに酸性雨・残留農薬の拡散・地球温暖化に伴う気候の変化・外来種による捕食など多くの要因が上げられている（松井, 1996；Blaustein et al., 2003; Pounds and Puschendorf, 2004）。

　あらゆる動物群において個体群サイズが毎年変動することは明らかであるが、その幅は個体群や動物群によって異なる。個体群変動の大きさは、その個体群の絶滅リスクにも影響すると考えられ、個体群の長期的な減少・増加傾向の検出のためにも必須の情報であることから注目を得るようになってきた（McArdle et al., 1990; Gaston and McArdle, 1994; Marsh, 2001）。両生類、特に陸生のサンショウウオ類は、トカゲ類とともに個体群変動が小さいとする報告もあるが（Gibbs, 2000）、日本産の両生類に関しては、そのような傾向を議論できるような長期的な個体群データは残念ながらまだ少ない。本報告では、著者が扱ってきたトウキョウサンショウウオ（Hynobius tokyoensis）やアカハライモリ（Cynops pyrrhogaster）などの両生類個体群の長期的モニタリングの結果を紹介すると同時に、既存の調査結果を幅広く収集し両生類の個体群サイズの年次変動の大きさについて概観した。また、著者らが行ってきたトウキョウサンショウウオ個体群の調査から得られた生活史・個体群パラメータを参考に、一般的な生活史をもつ両生類の個体群動態モデルを作成した（草野・川上, 1999; Akҫakaya et al., 2002）。動態モデルにより様々な条件下で長期間の個体群動態を予測し、その結果より個体群パラメータと個体群変動の大きさの関係を解析し、両生類の個体群変動の大きさに影響を与える個体群パラメータを探索した。

 

方法

繁殖個体群のモニタリングデータ

　繁殖個体群の調査は、東京都において2箇所で行ってきた。日の出町羽生の一つの繁殖池においては、1月下旬から4月に掛けてトウキョウサンショウウオの繁殖個体群を1975年より現在まで継続して毎年モニタリングを行ってきた（Kusano, 1980, 1982; 草野・井上, 2006; Kusano and Inoue, 2008）。繁殖メスの個体数を推定するため、繁殖シーズンにほぼ週1回程調査地を訪れ、産卵されている卵嚢の数を記録した。その際、卵嚢ごとに卵の発生段階にも注意して新しく産卵された卵嚢を特定した。繁殖メスは1対の卵嚢を産み出すので、総卵嚢数の半分がメスの個体数となる。1991年からは、八王子市南大沢のイモリ池でも調査を開始し、アカハライモリ・ヤマアカガエル（Rana ornativentris）の自然個体群とともに、移植個体群であるトウキョウサンショウウオとモリアオガエル（Zhangixalus arboreus）のモニタリングを行った（草野, 2014）。各種の繁殖シーズンにほぼ毎日池に通い、新しく産卵されている卵塊をカウントしその年の繁殖に参加したメスの個体数を推定した。アカハライモリに関しては、1メスが1個ずつ長期間に渡り産卵するため卵塊による調査は不可能なので、標識再捕調査により繁殖メスの個体数推定を行った（草野, 2014）。4月から11月までの活動シーズン中に原則として月一回の夜間センサスを行い、腹面の色彩パターン及びPITタグなどによる個体識別を行った（草野, 2008）。個体の捕獲歴データから、開放個体群モデルであるCormack-Jolly-Seber法を適用して個体数を推定した（Amstrup et al., 2005）。その他、日本産在来種だけでなく海外の両生類個体群のモニタリング結果の時系列データも文献より収集した。

 

年変動の大きさの指標

　動物の個体群は時間的に（及び空間的に）幅広い動態を示すことはよく知られているが、この動態を広く種間・個体群間で比較するためには適切な指標を用いる必要がある。個体数の年変動の大きさの指標としては、様々な指標が今までに提案され用いられているが、いずれも様々な問題点が指摘されている（McArdle et al., 1990; Gaston and McArdle, 1994）。その中で比較的よく使用されてきた未変換の個体数の変動係数（CV）を、今回変動の大きさの指標として用いた。ただし、様々な個体群間の比較には、多くの個体群が示す増加・減少などの長期的トレンドの調整が必須と考えられる。また、調査期間の長さも当然CVなどの変動の指標にも影響する（e.g., Marsh, 2001）。今回解析に用いた時系列データの長さは5〜42年に渡り、これら長さの違う時系列データを比較するためにGibbs（2000）の手法に従い指標を求めた。すなわち、直線回帰の残差を用いて長期トレンド傾向の調整を行い、かつ移動平均の手法に準じて、5年単位でCVを求めその平均を調整済み変動係数（detCV）として変動の大きさの指標とした。 今まで頻繁に使われていた上記の指標はいずれも平均値からの偏差に基づき変動の大きさを評価するもので、分布の正規性が問題となることがある。それに対して、近年提唱された指標である比例変動性（PV; Proportional Variability）は、毎年の個体数の対のすべての組み合わせについて、その比を求め平均することで変動の大きさを評価するので、分布の型には関わらず幅広く比較ができるノンパラメトリックな測度であるという（Heath, 2006; Heath and Borowski, 2013）。そこで今回の解析でも、変動の大きさの指標として前述したdetCVに加え、PVも参考のため示した（ただし、PVはトレンドや調査期間の長さに関して未調整）。

 

個体群動態モデルによる変動の大きさに影響を与える要因の解析

　個体群間の変動の大きさの違いの原因を探索するために、個体群動態モデルを作成しコンピュータ上で個体群変動のシミュレーションを試みた。個体群動態モデルは主に保全生態学の観点から、希少種の絶滅確率などを推定するための個体群存続可能性分析に利用されてきた（Reed et al., 1998; 鷲谷, 1999）。日の出町羽生におけるトウキョウサンショウウオの個体群研究によって得られた生活史・個体群パラメータを参考にし、一般的な生活史をもつ両生類の個体群動態モデルを構築した（Kusano, 1980, 1981, 1982; 草野・川上,1999; Kusano et al., 2006）。トウキョウサンショウウオなど多くの両生類の個体群は年齢構造を持つので、今回の解析ではレスリー行列を用いた齢構造モデルを採用し、個体数に影響する要因として人口学的確率性・環境の確率性・密度依存性・カタストロフィを導入した（Reed et al., 1998; 鷲谷, 1999; Akҫakaya et al., 2002）。

　個体群サイズはメス個体のみに注目して、繁殖メス個体数の年次変動を予測した。密度依存性としてはcontest型（Beverton-Holt型）の関数を用い、カタストロフィは数年に一度すべての年齢群の個体数が同時に半減するという過程を組み込んだ（Akҫakaya et al., 2002）。動態モデルには計10個の個体群パラメターが組み込んだ。なお、今回の動態モデルのプログラミング、シミュレーションの実行及び結果の統計解析には、統計言語Rバージョン3.5.1を利用した（R Core Team, 2018）。 McCarthy et al. (1995) の手法に準じて、動態モデルに組み込んだ10の個体群パラメータにランダムに数値を割り当て、その後の30年間の動態をシミュレーションした。それぞれのパラメータは、トウキョウサンショウウオの個体群パラメータの実際の観察値を参考に範囲を設定し、一様乱数を発生させて得た値を利用した）。一組のパラメータに対してシミュレーションを10回繰り返し、その度に変動の大きさの指標（detCVとPV）を求め、そのパラメータセットに対しての変動の大きさは10回の試行で得られた値の中央値を用いた。この作業を10,000回繰り返し、様々な個体群パラメータの組み合わせに対する変動の大きさを推測した。 シミュレーションによって得られた結果から、10個の個体群パラメータと30年間の年変動の大きさの関係を解析した。

　このような解析に良く用いられる一般化線形モデルでは、非線形の関係が想定されるデータの解析には必ずしもふさわしくない。そこで、柔軟性と予測精度の高さから幅広く予測モデルに使われている機械学習の一手法ランダムフォレスト・モデルを適用した（Breiman, 2001; James et al., 2018）。変動の大きさを目的変数にとり、10個の個体群パラメータ（表1参照）を説明変数としてRパッケージrandomForestの関数を利用して解析した（Liaw and Wiener, 2002）。モデルの予測精度は、モデルの作成に用いずに取り置いていたOOB（Out-of-Bag）データにより評価し、各個体群パラメータの変動の大きさへの影響は、重要度と部分従属プロット（Partial Dependence Plot）により評価した。重要度は、各変数をランダムに並び替え、それによる予測精度の減少の程度で評価する指標を用いた（Liaw and Wiener, 2002）。

 

　結果

両生類個体群の年変動の大きさ

　今回の解析では、在来種7種17個体群を含む総計16種35個体群の時系列データを用いて個体群年変動の大きさを評価できた（表1）。変動の大きさは、種・個体群によって大きく異なった。今回用いた2種類の指標（detCVとPV）では、トレンドや調査期間の長さの調整の有無が関係するのか、微妙に違いも認められたが、両指標間で強い正の相関が観察され（r = 0.791, P << 0.001）、どちらでもほぼ同じような傾向が見られた。

　表1を見ると、有尾類の種で変動が少ない傾向が見られた。変動係数（detCV）が0.3未満の個体群は、有尾類では 8例中 5例見られたのに対して、無尾類では 27例中 3例に過ぎなかった。特にトウキョウサンショウウオ（日の出町）とアカハライモリ（八王子）は、調査期間が30〜40年間と比較的長いにも関わらず、0.1台と極めて小さな値を示し、安定しているように見えた。これら個体群年変動の実際の様子を見るために、表1のうち調査期間が20年以上の時系列データを図1示した。また、Gibbs(2000)に報告されている両生類32種38個体群の調整済み変動係数のデータを加え、両生類の個体群変動の大きさを概観した（図2上図）。変動係数は0.05〜2.8を示し、平均値は0.633（SE = 0.054, N = 73）であった。図2には、参考のためGibbs(2000)の他の動物群（魚類・爬虫類・鳥類・哺乳類・蜘蛛型類・昆虫類）のデータもプールして示し、平均値はそれぞれ0.571、0.438、0.517、0.513、0.643、0.834であった。動物群間では、昆虫類のみ他より有意に高い傾向が見られたが（t検定による対比較; Holm法による多重検定調整済みP < 0.05）、昆虫類以外の動物群間では有意な差は検出できなかった。両生類以外の動物群でも、変動係数の分布はほぼ両生類と同じ傾向が見られ（図3下図）、トウキョウサンショウウオとアカハライモリの示した個体群変動の小ささは、動物群を広く見渡してみても際立つことが分かった。

 

図1.　両生類の個体群サイズの年次変動の実例。それぞれのデータの引用元は表1 参照。

 

変動の大きさに影響する要因

　ランダムフォレスト・モデルを適用することにより、シミュレーションで得られた変動の大きさの値へ良く適合した予測モデルを得ることができ、目的変数（detCVとPV）の分散のそれぞれ83％と89％をモデルにより説明できた。その結果、変動の大きさへの影響力の強さは、個体群パラメータにより異なり、変態後年生残率が他に比べて極めて強く影響していることが分かった（図3）．逆に、カタストロフィの間隔・変態後年生残率のCV・繁殖率のCVなどはほとんど関与していないことが分かった。 影響が示唆された7個のパラメータへの変動の大きさの反応パターンを、部分従属プロットにより示した（図4）。変態後年生残率以外の6個のパラメータでは、パラメータの値が増加すると変動の大きさの予測値はすべて単調な増加・減少傾向パターンを示したが、変態後年生残率では0.7〜0.8までは予測値が急激な減少を示し、それ以降は一定レベルを維持し飽和するパターンを示した。今回の解析では、変動の大きさにどの指標（detCVとPV）を使っても、初期個体数を除き、ほぼ同じ結果を得ることができた（図4）。

図2.　両生類個体群の変動の大きさの頻度分布。Gibbs（2000）のデータに今回収集した文献データの結果も加えて集計した。比較のために，両生類以外の動物群（魚類・爬虫類・鳥類・哺乳類・昆虫類・蜘蛛類）も集計して示した。

図3.　ランダムフォレスト・モデルによる変動の大きさに影響を与える要因の解析結果。変動の大きさの指標として，調整済み変動係数（detCV）と比例変動性（PV）を用いた．詳細はテキスト参照のこと。

図4.　変動の大きさの各個体群パラメータへの反応曲線（部分従属プロット）。変動の大きさとして，調整済み変動係数（detCV）と比例変動性（PV）の結果をそれぞれ示した。

 

　 考察

　今回の結果は、一般的に動物個体群の年変動の程度は、種や個体群により大きく変異することを示した先行研究の結果を支持するものだった（Gibbs, 2000; Marsh, 2001; 表2）。両生類の中で最も変動が大きかったのは、スキアシガエル（南カロライナ）で変動係数が2.78、最も小さかったのは、コオロギガエル（フロリダ）で変動係数が0.05だった（Gibbs, 2000）。モリアオガエル（大峰沼）、トウキョウサンショウウオ（日の出町）、アカハライモリ（八王子）と在来の3個体群がこれに続きいずれも0.1台の小さな値を示した（表1）。トウキョウサンショウウオは移植個体群（八王子）でも、移植からその後の長期にわたる定着過程でも0.2台と小さな値を示している（表2; 草野, 2014）。一方、モリアオガエルの他の3個体群は0.36〜0.46と必ずしも種として変動が小さいとは言えないかもしれない。

　両生類個体群を想定した動態モデルによるシミュレーションの結果から、個体群の年変動の大きさには様々な個体群パラメータが関与していることが推察された（図3・4）。特に、変態後の年生残率が大きな影響を与え、生残率が低いほど個体群は大きく変動し、生残率が高いと安定する傾向があることが分かった。トウキョウサンショウウオ（日の出町）に関しては、標識再捕の結果より変態後の年生残率は0.7と推定され（Kusano, 1982）、また年齢分布からもほぼ同じ推定値が得られている（草野他, 2014;未発表データ）。同様に変動の小さかったアカハライモリでも、標識再捕の結果より0.7〜0.9と極めて高い変態後の年生残率の推定値が得られている（草野, 2014;未発表データ）。逆に、CVの値が0.642と高かったナガレタゴガエルでは、繁殖メスの年生残率は繁殖個体の年齢分布から（Kusano et al., 1995）0.3と低目の値が推定されている（未発表データ）。一方、CVの値は低くはないがそれほど高くもない0.3〜0.4台の値を示しているニホンアカガエルでは、東広島の個体群でメス成体の年生残率が0.48と推定されている（Marunouchi et al., 2002）。これら個体群の年生残率に関する情報は、動態モデルのシミュレーションによる解析結果の妥当性を支持するものかもしれない。もしそれが正しいとすれば、CVの値が0.9〜1.0を越えたヤマアカガエルの2個体群などは（表2）、成体の年生残率がかなり低い可能性が示唆される。いずれにしても、解析結果の妥当性の検証は今後の課題である。

　動態モデルによるシミュレーションでは、変動の大きさだけでなく30年間の絶滅率も評価できる。この両者間には有意な強い相関が見られ（r = 0.794, N = 10,000, P << 0.001）、予想されるとおり変動の大きな個体群は絶滅するリスクも高い傾向がある。つまり、変態後の年生残率の低い個体群は、変動が大きいだけでなく絶滅リスクも高いということになる。一般に止水で繁殖する両生類では、変態までの死亡率が極めて高い傾向が強いが（e.g., Kusano, 1981）、卵・幼生期よりもその後の変態後の生残過程が、個体群の絶滅にも強く関わるかもしれない（図3）。したがって、個体群の保全を考える際繁殖池の環境は当然だが、それ以上に幼体・成体が生息する陸上の環境も重視する必要があるだろう。特に、繁殖個体の死亡・損失は個体群の存続に重大なインパクトを与えるかもしれない。近年問題となってきた外来種による捕食や販売目的の繁殖個体の採集は、両生類個体群の保全にとって対処すべき最重要課題の一つとなってきている。野外での長期間の個体群調査と、その結果の詳細な解析の必要性は今後ますます高まるに違いない。今後のこの分野の研究の進展に期待したい。   　

 

引用文献

　

この小文は、下記の報文を基に作成したものです。詳しくは下記の元論文をご覧ください。別刷りPDFが必要な方はここに請求してください。

　
草野　保．2019．両生類の個体群変動の大きさに影響を与える要因．p.96–106．日本両生類研究会 20 周年記念誌編集委員会（編）両生類に魅せられて：カエルとサンショウウオの長期研究と最新の研究成果．日本両生類研究会，新潟．

WB01343_.gif (599 バイト)戻る