Sampling
(i)Genome-wide DNA data extracts from teeth of skull (MCE1356, Fig.) を解析。 2b, Supplementary Fig. S1）、および西グリーンランドのディスコ湾で採取したベルーガ8頭とイッカク8頭の組織サンプルから抽出したゲノムワイドDNAデータと、（ii）MCE1356と西グリーンランドのベルーガ18頭とイッカク18頭からなる参照パネルの骨コラーゲンの安定炭素・窒素同位体組成（図1c）の解析を行った。 組織サンプル（皮膚）は96%エタノールで保存されました。 これらのサンプルは、グリーンランド政府からのイヌイットの生物学的サンプリングに関する一般許可に基づき、グリーンランド天然資源研究所の科学者が収集し、CITES許可証IM 0401-897/04, IM 0721-199/08, IM 0330-819/09, 116.2006 によりデンマークに輸出されたものである。 サンプル情報は補足表S2に詳述。
DNA extraction, library preparation and sequencing
Tissue samples
DNA from tissue samplesはQiagen Blood and Tissue Kitを用いてメーカーのプロトコルに従い抽出された。 Covaris M220 Focused-ultrasonicator を用いてDNAを断片化し、約350-550塩基対（bp）の断片長を作製した。 イルミナNeoPrepを使用し、NeoPrep Library Prep System Guideに従って、デフォルト設定を適用して、断片化したDNA抽出物からライブラリーを構築した。 PCR増幅、定量、正規化はすべてNeoPrep Library Prep Systemによって行われた。 ライブラリーは、Agilent 2100 Bioanalyzerでサイズ分布をスクリーニングし、等モル比でプールしてから、80 bp SEテクノロジーのIllumina HiSeq 2500でシーケンスした。
歯/骨サンプル
組織サンプルと異なり、古い骨や歯のサンプルは比較的DNA濃度が低く、異なる抽出およびライブラリー構築プロトコルが必要となる。 標本MCE1356の5本の歯と1つの骨片から約0.5gの骨粉を、手持ちのドレメルを使って穴を開けた。 コペンハーゲン大学デンマーク自然史博物館の古代 DNA クリーンラボで、11 に記載した抽出バッファーに 30 分の消化前段階12 を加えて、歯や骨の粉末から DNA を抽出した。 Zymo-Spin V カラム (Zymo Research) の代わりに、Amicon Ultra 30 kDa Centrifugal Filter Units を用いて抽出バッファーを濃縮し、Qiagen Minelute チューブを用いてさらに濃縮と洗浄を行った。 精製された抽出物は、次に13によって記述されたプロトコルに従ってイルミナライブラリーに構築された。 qPCRを使用して、ライブラリ構築が成功したことを確認し、配列決定するライブラリを選択し、過増幅を引き起こすことなく各ライブラリを十分に増幅するために必要な適切なPCRサイクル数を計算した。 合計で4つのライブラリーがユニークな6bpインデックスで増幅され、250bp SE配列決定を用いてIllumina MiSeqプラットフォームで内因性コンテンツをスクリーニングした。 ベストライブラリーは、次に80bp SE技術を使用してIllumina HiSeq 2500プラットフォームで再シーケンスされた。
バイオインフォマティクス解析
すべてのマッピングおよびDNA損傷分析は、Paleomixパイプライン1.2.1214内で行われた。 リードはAdapterRemoval 2.2.015で、最小リード長を25bpに設定した以外はデフォルトの設定でトリミングした。 リードはFastQCで検査し、BWA16でBacktrackアルゴリズムを適用してアライメントし、開始種長を無効にした。 複数の位置にマッピングされたリードやマッピング品質スコア（BWAによるMAPQスコア）が30未満のリードは、SAMtools17を使用して削除しました。 配列の重複はPicardのMarkDuplicates (available from: http://broadinstitute.github.io/picard) を用いて除去し、最終アライメントはGATK18を用いてインデルの周囲で再整列させた。 検体MCE1356のシトシンからウラシルへの脱アミノ化は、mapDamage v2.0.619の出力を用いて評価した。 MapDamageの結果、配列に明確な脱アミノ化のシグナルは見られなかった（補足図S3）。
ミトコンドリア解析
MCE1356の母系を決定するため、DNA配列決定リードをベルーガ（GenBank accession: KY444734）とイッカク（GenBank accession: NC_005279）のミトコンドリア参照ゲノムと対応させて平均範囲を比較検討した。 参照パネルを構成する8つのベルーガと8つのイッカク試料からのリードは、それぞれのミトコンドリア参照ゲノムにマッピングされた。 BEDtools20、SAMtools17、GATK18を用いて、5リード以上でカバーされる領域を持つMCE1356と16の参照パネルサンプルのミトコンドリアコンセンサス配列を構築した。 ClustalW21を用い、16の参照パネルサンプルとベルーガのミトコンドリア参照ゲノムにマッピングされたMCE1356のバージョン、またはイッカクのミトコンドリア参照ゲノムにマッピングされたMCE1356のバージョンを含む、デフォルト設定を適用した2つの配列アライメントを作成しました。 この2つのアラインメントを用いて、Popart 1.723 (available from: http://popart.otago.ac.nz) を用いて、インデルのある部位や欠損データを除いて、中央結合型ハプロタイプネットワーク22を構築した。 その後、マッピング後のカバレッジと2つのハプロタイプネットワークの両方を用いて、MCE1356の母系種を決定した。
Nuclear DNA analysis
すべてのサンプルから得られたDNAシーケンスリードをシャチ（Orcinus orca）参照ゲノム（GCA_000331955.2）にマッピングした。 最近、高品質のシロイルカのゲノムが発表されたが24、2つの親となりうる種のうちの1つにマッピングすると、解析に偏りが生じる可能性がある。 また、シャチはベルーガやイッカクと比較的近縁でありながら（分岐時間12 MYA）25、我々の解析を複雑にする内生殖のリスクを最小化するのに十分な距離であることから、シャチゲノムにリードをマップした。 ANGSDに実装されているSAMtools17を使用して遺伝子型尤度を推定し、27に記載されているように最尤法を使って遺伝子型尤度から直接メジャーアレルとマイナーアレルを推測しました。 全個体のデータを合わせて可視化することで、平均読み取り深度を4.14倍と推定し、読み取り深度が高い部位を特定することができた。 このような部位は、ゲノムのパラログや繰り返し領域に由来している可能性が高いことが分かりました。 データセットを視覚的に検査し、さらにフィルタリングして、読み取り深度が9を超えるすべてのサイト（サイトの6.9%）を除外しました。 ANGSDでは、SNPはその対立遺伝子頻度に基づいて呼び出された。 マイナーアレル頻度（MAF）は遺伝子型尤度から推定し、尤度比検定を用いてMAFが0と異なるかどうかを判断した。 尤度比検定のp値が<1e-4であれば、その部位は多型とみなされ、データセットに保持された。 このSNPフィルターを適用すると、4本以下のリードでカバーされる部位は、これほど低いp値を得ることができないので、保持されないことになる。 これらのフィルターを、17個体すべてを含むデータセットと、MCE1356を含まないデータセットに適用した。
MCE1356がベルーガとイッカクの雑種かどうかを判断するために、さらに17個体のデータセットにフィルターをかけ、MCE1356のリードを持たない部位を除いた。
我々の目的はMCE1356で見つかった対立遺伝子を参照パネルの対立遺伝子と比較することだったので、異なる参照パネルの最小ユニークリード深度と対立遺伝子頻度を与えてMCE1356で見つかった対立遺伝子を間違った親種に割り当てる確率を推定した。 ユニークリード深さは、特定の部位をカバーするリードの数として定義され、すべてのリードはユニークな個体から得られたものである。 この確率Pは、式1のように計算された。
ここでf（ps）は親種における対立遺伝子周波数、url（psp）は参照パネルにおける種固有の固有リード深さである。
最高の確率f(ps-max)を与える親種の対立遺伝子頻度は、式2のように記述することができる。
式1に式2を挿入して、式1, MCE1356で見つかった対立遺伝子を誤った親種に割り当てる最大確率P(max)は、式3のように計算された。
結果、ユニークリード深さが2, 各親種で3つ、4つの場合、MCE1356で見つかった対立遺伝子を間違った親種に割り当てる最大確率は0.148、0.105、0.082であった。 これらの最大値は、すべての部位が親種内で変動し、すべての部位のMAFがちょうど（1 – (urd/urd + 1)）である場合にのみ得られるので、誤差率と解釈すべきではない。 さらに、ベルーガとイッカクのMAF分布が類似していると仮定すると、MCE1356で見つかった対立遺伝子の誤った割り当てが両種に等しく影響し、したがって、交雑の推論に最小限の影響を与えることになるであろう。 式2、3においてユニークリード深度を用いる利点は、MCE1356で見つかった対立遺伝子が誤った親種に割り当てられる確率を推定する際に、個体数とリード数を組み合わせることができる点である。 その後、リード深度分布、MAF（メジャーアレルとマイナーアレルを固定した場合）、各サイトにデータを持つ個体数を親種ごとに別々に計算した。
MCE1356以外に、最小ユニークリード深度が2、3、4の親種パネルを含む3つのデータセットを構築した。 3つのデータセットで保持されたサイト数は、それぞれ61,105、11,764、360であった。 部位数を最大にすることと、MCE1356で見つかった対立遺伝子の誤った割り当てを最小にすることの妥協点として、MCE1356で1リード、各親生物種で最小ユニークリード深さ3のデータセットを使用することにした。 そのデータセットには11,764部位が含まれ、以降の解析に使用した。
17個体のデータセットについて、(i) ベルーガとイッカクの種パネルで異なる対立遺伝子に固定された部位数、 (ii) ベルーガで多型であるがイッカクで多型ではない、 (iii) ベルーガではなくイッカクで多型、 (iv) ベルーガとイッカクの両方で多型などの要約統計処理を実施した。 2つの親種のパネル間で固定されていると推定される部位は、2つの親種間で高度に分化している、すなわち、アリル頻度に大きな差があるマーカーが濃縮されていることになる。 これらのマーカーは、必ずしも2つの親種間で固定された差ではないものの、MCE1356の家系に対して高い情報量を持つ。
MCE1356を含まないデータセットからの遺伝子型尤度を用いて、参照パネル中のベルーガとイッカクが、それ自体最近混血した個体ではないことを検証した。 NgsAdmix28を用いて、2つの集団（K = 2）を指定しながら、個々の混血係数を推定した。 100回実行し、混血係数の平均値と標準偏差をその後の解釈に使用した。
我々は、fastNGSadmix29を用い、1000ブートストラップを適用してMCE1356の混血比率を解析した。fastNGSadmixは、参照集団/種の対立遺伝子頻度と単一個体の遺伝子型尤度を用いて、その混血比率を推定している。 このソフトウェアは、0.00015×29という低いカバレッジのNGSデータで堅牢であることが証明されており、したがって我々の研究に最適であった。
我々はさらに、ベルーガパネルとイッカクパネルで異なる対立遺伝子について固定した部位（9178部位）を調べ、MCE1356のベルーガ特異的対立遺伝子とイッカク特異的対立遺伝子に一致するリードの観測割合を異なる混成シナリオの下で期待される割合と比較し、MCE1356の交雑状況を推定した。 7つの異なるハイブリダイゼーションシナリオが、どの程度、観測されたデータと一致しているかを判断するために、ピアソンのカイ二乗適合度統計量を計算した。 この検定統計量は式4のように計算される：
ここでOiとEiはそれぞれ親種i由来のアレルの観察数と期待数である。 選択したシナリオが観測データによく対応するという帰無仮説の下では、検定統計量Tは中心χ2分布に従う。 したがって、観測データと最もよく対応するシナリオが最も低い検定統計量を導くことになる。
7つの異なる混成シナリオをさらに調べるために、親集団の異なる対立遺伝子に固定された部位で観測される対立遺伝子の尤度を計算した。 具体的には、親種からの対立遺伝子を二項モデルで継承し、さらにマーカー間の独立性を仮定して、雑種MCE1356で観察される対立遺伝子の尤度を計算した。 しかし、独立性の仮定を破ってもなお、不偏の推定値が得られることに注意したい。 雑種MCE1356がベルーガ祖先の割合bとイッカク祖先の割合（1-b）からなると仮定すると、bの尤度は、部位iでベルーガ対立遺伝子と一致するリードの数nibとイッカク祖先と一致するninが与えられたとき、式5のように独立したk部位での尤度の積として書くことができる。
bの対数尤度をその全範囲にわたって計算した。 から1まで、そして7つの異なる混成シナリオで対数尤度を比較した。
Sex determination
MCE1356 とシロイルカとイッカク参照パネルの個体の性別を決定するために、X染色体と常染色体のカバーレートを調査しました。 これは、X染色体に由来すると思われるスキャフォールドのカバレッジと、常染色体に由来すると思われるスキャフォールドのカバレッジを比較することによって行われたものである。 シャチゲノムと牛（Bos taurus）のX染色体（CM008168.2）をSatsumaSynteny30でアライメントし、どのスキャフォールドがX染色体に由来するかを決定した。 さらに、X染色体とY染色体の一部の領域の相同性に起因するバイアスを除去するために、得られた推定X染色体スキャフォールドをヒトY染色体（NC_000024.10）にさらにアライメントし、これらのスキャフォールドを以降の解析から除外した。 次に、SAMtools17の深さ関数を用いて、残りのXスキャフォールドとX染色体に整列しない4つの大きなスキャフォールド（すなわち常染色体スキャフォールド）全体のカバレッジを計算した。 ゲノム全体のカバレッジのばらつきを補正するために、10,000,000サイトをランダムに選択し、これらのサイトの平均カバレッジを計算し、このステップを100回繰り返した。 各個体について、5%および95%信頼区間、ならびに第1および第3四分位を算出し、その後の解釈に用いた。
安定炭素・窒素同位体分析
MCE1356、ベルーガ18頭、イッカク18頭の頭蓋骨から粉末状の骨試料（〜100 mg）を10 ml 2:1 クロロフォルム・メタノール（v/v）で超音波処理下に1時間処理して、脂質を除去した。 溶媒を除去した後、常圧で24時間乾燥させた。その後、試料を10 mlの0.5 M HClで4時間、オービタルシェーカーで撹拌しながら脱塩した。 脱塩後、試料をI型水にて中性まで洗浄し、10-3M HCl中75℃で36時間加熱し、コラーゲンを可溶化させた。 その後、水溶性コラーゲンを凍結乾燥した。 次に、このコラーゲン抽出物を10mlの10：5：4クロロホルム／メタノール／水（v／v／v）で1時間超音波処理し、残留する脂質を除去した。 遠心分離後、クロロホルム/メタノール層を除去し、60℃で24時間、水層からメタノールを蒸発させた。サンプルを再び凍結乾燥し、元素および同位体分析用に錫のカプセルに秤量した。 炭素と窒素の安定同位体組成と元素組成は、カナダのトレント大学でNu Horizon (Nu Instruments, UK) 連続フロー同位体比質量分析計を用いて測定した。 炭素と窒素の安定同位体組成は、USGS40 と USGS41a を用いて VPDB と AIR スケールに照らして校正された。 標準不確かさはδ13Cで±0.17‰、δ15N31で±0.22‰と決定され、追加の分析の詳細は補足資料S4に記載されている。 ベルーガとイッカクの同位体組成の違いの統計的有意性は、対応のないt検定を用いて評価した。