Convergent genomic signatures of domestication in sheep and goats

Posted by admin Articles

Sampling

Iran (IROA, IRCH group) と Morocco (MOOA, MOCH group) で家羊 (O. aries) と山羊 (C. hircus) を1グループ20匹ずつサンプルしました（補図6）。 これらの試料は，欧州連合指令 86/609/EEC の倫理規定に基づき，Nextgen European プロジェクト（助成金契約番号 244356）の枠内で，モロッコ北部では 2008 年 1 月から 2012 年 3 月の間に，イラン北西部では 2011 年 8 月から 2012 年 7 月の間に採取されたものである． 4804>

家畜化ゆりかご内のイラン北西部で野生種Asiatic mouflon (O. orientalis) とBezoar Ibex (C. aegagrus) を採取した21,22。 13 個の Asiatic mouflons と 18 個の Bezoar ibex の組織（それぞれ IROO と IRCA グループ、補足図 6）は、飼育下または最近狩猟した動物から、またイラン環境局で入手できる冷凍試料から採取したものである。 この個体ベースのサンプリング手法は、局所的な影響（例：, 4804>

追加データ

さらに、羊と山羊の世界品種パネル（それぞれwpOAとwpCH）を組み立てた。wpOAは、International Sheep Genome Consortiumが提供する20種類の世界品種を表す12倍範囲の全ゲノム再シーケンス（WGS）サンプルで構成されていた。 フランスアルペン種2個体、フランスサーネン種2個体（INRAによる配列決定）、イタリアサーネン種5個体（Parco Tecnologico Padanoによる提供）、オーストラリア5個体（例, 4804>

WGSデータの作成

ゲノムDNAは、Macherey Nagel NucleoSpin 96 Tissue Kitを使用して、メーカーのプロトコルに適合させ、すべての組織サンプルから抽出することに成功した。 組織サンプリングは、1サンプルあたり25mgの断片を得るために、MNスクエアウェルブロックで行った。 3.5個のMNスクエア-96ブロックを用意し、Tecan Freedom Evo Liquid handlerを使用して、製造元のプロトコルにしたがって抽出を実施した。 180 µl の T1 Buffer と 25 µl の proteinase K を用いて、56 ℃で一晩サンプルをホモジナイズする予備溶解ステップを実施した。 結合条件を調整するために、200 µl のBQ1 バッファーを加え、サンプルプレートを70 °Cで1時間インキュベートし、その後、200 µl の100%エタノールを添加した。 ライセートをNucleospin Tissue binding plateに移し、真空（-0.2 bar、5分）をかけてフロースルーを除去した。 BWおよびB5バッファでそれぞれ3回洗浄し、再び真空を適用してフロースルーを廃棄した。 ゲノムDNAの溶出に先立ち、Nucleospin Tissue binding plateシリカメンブレンを少なくとも-0.6 barの真空下で10分間乾燥させた。 溶出ステップは、96-PCRウェルにあらかじめ温めた100 µlのBEバッファー（70℃）と3700 rcfで5分間の遠心分離ステップで行った。 ゲノムDNAは凍結融解を避けるため4℃で保存し、ピコグリーン法とナノドロップを用いて濃度（ng/μlとして）検査を行った。

サンプルごとにCovaris® E210装置を用いて150～700 bpの範囲でせん断したゲノムDNA 500 ngから全ゲノムを再スクリーニングし、半自動プロトコルによりイルミナ®ライブラリー作製に使用した。 End repair、A-tailing、Illumina®互換アダプター（BioScientific）のライゲーションは、SPRIWorks Library Preparation SystemとSPRI TE instrument（Beckmann Coulter）を用いて、メーカーのプロトコルに沿って実施しました。 300-600 bpのサイズセレクションを行い、ほとんどの断片を回収した。 Platinum Pfx Taq Polymerase Kit (Life® Technologies) とIllumina®アダプター特異的プライマーを用いた12サイクルのPCRにより、DNA断片が増幅された。 ライブラリーは0.8x AMPure XP beads (Beckmann Coulter) で精製し、Agilent 2100 Bioanalyzer (Agilent® Technologies) とqPCR定量で解析した。 ライブラリーは、Illumina® HiSeq2000のペアエンドフローセルで100塩基長のリードケミストリーを用いてシーケンスした。

OvisのIlluminaペアエンドリードは、BWA-MEM48を用いて羊参照ゲノム（OAR v3.1, GenBank assembly GCA_000298735.146) に、Capraについてはヤギ参照ゲノム（CHIR v1.0, GenBank assembly GCA_000317765.147) にマッピングされた。 各個体について作成されたBAMファイルはPicard SortSamでソートされ、逐次Picard MarkDuplicates（http://picard.sourceforge.net）、GATK RealignerTargetCreatorおよびGATK IndelRealigner49、SAMtools calmd50を用いて改良した。

Variant discoveryには3種類のアルゴリズムを用いて実施された。 4884>

変異体の発見は、SAMtools mpileup50、GATK UnifiedGenotyper51、Freebayes52の3つの異なるアルゴリズムを使用して行われました。 変異部位は、呼び出しアルゴリズムのマルチサンプルモードを使用して、6つのグループのそれぞれについて独立して同定された。 (i) MOOAから162サンプル、(ii) IROAから20サンプル、(iii) IROOから14サンプル、(iv) MOCHから162サンプル、 (v) IRCHから20サンプル、 (vi) IRCAから19サンプル。 いくつかのグループについては、NextGenプロジェクトの一環として、より多くの個体のWGSが利用可能であった（上記参照）。 本研究で使用したサンプルは、可能な限り20個体のバランスのとれたグループを得るために選択された。 IRCAおよびIROOグループについては、後の段階で追加サンプルが入手可能となり、下流の解析のために追加した。 アラインメントとコーリングの質が低い動物は、最終データセット（補足データ5）を得るために削除された。

各グループ内で、バリアントサイトの質のフィルタリングが2ラウンド行われた。 Filtering Stage 1では、3つのアルゴリズムからのコールをマージし、最も信頼度の低いコールをフィルタリングした。 少なくとも2つの異なるコーリングアルゴリズムで、phred variant quality >30でコールされた場合、variant siteはパスしました。 ある部位の代替アレルは、コーリングアルゴリズムのいずれかによってコールされ、遺伝子型カウントが>0であれば合格。 フィルタリング第2段階は、GATKによるVariant Quality Score Recalibrationを使用しました。 まず、(i)その部位が3つのバリアントコールャーのすべてによって、phred-scaled variant quality >100でコールされ、(ii)その部位がバイアリルで、(iii) genotype phred-scaled quality >30でサンプルだけをカウントしながらマイナーアリルカウントが少なくとも3であるグループ内の最も信頼度の高いバリアントサイトのトレーニングセットを生成しました。 トレーニングセットは、UnifiedGenotyperから以下のバリアントアノテーションを使用して、ツールGATK VariantRecalibratorを使用してガウスモデルを構築するために使用された。 QD, HaplotypeScore, MQRankSum, ReadPosRankSum, FS, DP, InbreedingCoefficient. データセット全体にガウシアンモデルを適用し、VQSLOD（真のバリアントであることの対数オッズ比）を生成。 VQSLOD <cutoff valueの場合、サイトはフィルタリングされた。 カットオフ値は、各グループで以下のように設定した。 最小VQSLOD = {トレーニングセットのバリアントのVQSLODの中央値}-3 × {トレーニングセットのバリアントのVQSLODの中央値絶対偏差}. OvisとCapraの6つのグループ（イランとモロッコの家畜、各属の野生種）のSNPsコールセットを作成し、transition/transversion SNP比から、選択性と感度のバランスが最も良いカットオフ基準を選択したことが示唆された。 本研究で行った解析ではグループ間の比較が必要であったため、どのグループの動物でも一貫したSNP部位で遺伝子型コールセットを作成した。 各属について、その3つのグループからバリアントコール部位をマージし、バイアリル位置のみを欠損データなく保持した。 GATK UnifiedGenotyperでGENOTYPE_GIVEN_ALLELESオプションを使用して、対象となるすべての個体の各バイアリルSNP部位で遺伝子型を再コールした。 この段階で、ヒツジとヤギの世界品種パネルに属する動物（wpOAとwpCH）と、この段階で入手可能となった追加の野生サンプル（O. orientalis 4匹とC. aegagrus 4匹）に、個体リストを拡張した。 遺伝子型はBeagle 453によってグループ内で改良・段階化し、遺伝子型確率が0.95以下のところをフィルターで除外した。 最後に、本研究で用いた個体の異なるサブセット間で単型である部位をフィルタリングした（下記参照）。

OvisとCapra間で検出された選択のシグナルを比較するために、2つの参照ゲノム間でクロスアライメントを行った。 まず、Ensembl release 69 code base54のpairwise alignment pipelineを用いて、ヒツジ（OARv3.1）とヤギ（CHIR1.0）の参照ゲノムを整列させました。 このパイプラインは、DNAレベルでのアライメントにLastZ55を使用し、その後、両ゲノム内の位置に応じてアライメントされたブロックを連結する後処理を行います。 LastZのペアワイズアライメントパイプラインは、Ensemblがサポートするすべての種についてルーチンに実行されていますが、ヤギはまだEnsemblに含まれていません。 どちらかの種に偏ることを避けるため、2種類の種間アラインメントを作成した。 一方はヒツジを参照ゲノム、ヤギを非参照とし、もう一方はヤギを参照ゲノム、ヒツジを非参照としました。 この違いは、参照種のゲノム領域は非参照種の単一遺伝子座に一意に対応することを強いられるのに対し、非参照種のゲノム領域は参照種の複数の場所に対応することを許されることである。 基準ゲノムの染色体セグメントについて、非基準ゲノム上の座標を求めた。 最後に、一方の属で発見されたSNPについては、他方の属の参照ゲノムとの全ゲノムアライメントを用いて対応する位置を特定した（補足表6）。

遺伝子構造

グループ内の遺伝的多様性を表すために、VCFtools56を用いて、Ovisについて73個体の遺伝的変異要約統計量を算出した（すなわち。 13 IROO、20 IROA、20 MOOA、20 wpOA）、カプラについては72個体（すなわち、18 IRCA、20 IRCH、20 MOCH、14 wpCH）についての遺伝的変異の要約統計量を算出した。 測定した統計量は、各属および各グループ内の個体全体における多型変異の総数（S）、各グループ内の平均塩基多様性（π）、各個体の近親交配係数（F）であった。 各属内において、野生グループと各家畜グループとの差は、個体の近交係数および遺伝的負荷値については片側t検定、部位ごとのヌクレオチド多様性については両側Mann-Whitney検定を用いて検定した

各属内の4グループ間の全体的な分岐（すなわち, VCFtools56に実装されているWeir and Cockerham57に従った二重鎖SNPと平均加重ペアワイズFstを用いて、各属内の4グループ（すなわち野生、イランとモロッコの家畜、世界パネル）間の全体的な分岐を見積もった。 グループ間の遺伝的構造は、VCFtoolsを用いて連鎖不平衡 (r²) が0.2より大きいSNPを除去するためにデータセットを刈り込んだ後、クラスタリング法sNMF26で評価した。 連鎖不平衡 (r²) は、50SNPのスライディング・ウィンドウ内のSNPのペア間で計算され、r²が0.2より大きい場合、ペアごとに1SNPがランダムに削除された。 各sNMF分析について、同じクラスタ数（K）の5回の実行が、Kの値を1 から10として実行された。 しかし，個体がクラスタ間でどのように分割されたかを評価するために，異なるKの数の代替パーティションも探索された。

共有祖先と混血を分離するために，sNMFに使用した刈り込みデータセットを用いて，集団分割とその後の混血イベントを共同で推定するためにTreeMix27を実行した。 最尤推定をより精緻にするために、-globalオプションをつけてTreeMixを実行した。 また、野生個体と家畜個体との分岐を考慮し、TreeMixの木を根付かせた。 ジャックナイフのブロックサイズは-k 500 SNPsで、これはおよそ150 kbに相当し、ヒツジとヤギの両方で見られるLDの平均ブロックを超えている。 移動がない場合の最尤樹を作成し、次に移動イベントを追加して、モデルによって説明される分散の増分変化と個体間の残差値を調べた。 目的は、野生個体と家畜個体の間に高い残差値や移動のエッジが存在する可能性を検出することであった。 TreeMixによって同定された混血の可能性のあるベクトルの統計的関連性をさらに調べるために（補足表3）、ADMIXTOOLS suite58のqp3Popプログラムを用いて、遺伝的導入の公式検定として3集団検定f328を各グループの組み合わせについて算出した。 カプラについては、wpCH群をオーストラリア系品種、フランス系品種、イタリア系品種に分けた。

Demographic inference

各属について、MSMC2ソフトウェア25に実装されているMSMCモデルを用いて祖先の人口動態推論を行った。 MSMCはpairwise sequentially Markovian coalescent59に基づくが、入力として段階的ゲノム配列データのハプロタイプを用いる。 各解析には1グループから2個体、つまり4つのハプロタイプを使用した。 各解析は、別のランダムな2個体セット、すなわちグループごとの解析の複製に対して繰り返された。 入出力ファイルは、MSMCソフトウェアに付属するpythonスクリプト（https://github.com/stschiff/msmc-toolsにある）を用いて生成・解析した。 解析パラメータは、突然変異率を2.5×10-8に設定し、世代数を2年に設定した以外は、デフォルトのままである。 これらのパラメータを変化させ（突然変異率2.5×10-8と1.0×10-8、世代長2年と4年の組み合わせ）、家畜化期間の概算を行った（補足図2参照）。 まず、Libradoらの方法に従って、各個体の遺伝的負荷を、すべてのタンパク質をコードするゲノム位置にわたる劣化的な適性効果の合計として計算することにより、60。 簡単に言えば、進化的制約の代理として、46ウェイ哺乳類アライメント（http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP46way/placentalMammals/）のPhyloPスコアを使用したのである。 このアラインメントから、機能的制約の下で進化しているタンパク質コード部位を同定した（phyloP score ≥1.5）。 次に、OvisまたはCapraの各ゲノムについて、これらの部位が変異しているかどうかを調べた。 もし変異があれば、すべての変異部位のフィロPスコアを合計し、制約の多い部位の変異が全体の負荷推定値に比例して多く寄与するようにした。 このようにして、各ヒツジ／ヤギゲノムの負荷推定値が得られた。 最後に、部位ごとの平均負荷を得るために、解析した位置の総数で割った。 なお、ヘテロ接合部位における変異の優性係数のモデル化（例：劣性、中間、優性）を避けるために、ホモ接合部位を条件としたことは注目に値する。 第二に、野生と家畜のオヴィスグループにおける遺伝的な劇症負荷を遺伝子ごとに比較し、代替仮説を家畜が野生近縁種よりも負荷が大きいとして、ウィルコクソン検定を行った。p値は多重検定61で補正し、調整p値<0.05の閾値を適用した。 遺伝的負荷の有意な増加を示す遺伝子セットについて、WebGestaltを用いた遺伝子オントロジー・エンリッチメント解析を行った62,63。 参照ゲノムは遺伝子の注釈が乏しいため、我々の種とヒトやマウスとの間のシングルコピー・オーソログに依存した。 X 染色体由来の遺伝子はバックグラウンドセットから除外した。

選択シグナルの検出

家畜化に関連する選択のシグナルの検出には、3つのグループ（すなわち、イラン、モロッコの家畜グループ、各属の野生グループ）の少なくとも1つにおいてマイナーアレル頻度が0.10より大きいすべての二重対立SNPsを使用した。 家畜化過程に関連する選択のサインがすべての家畜に存在すると予想されるため、以下の一般的な戦略を採用した：各属について、野生グループと伝統的に管理されている家畜グループ（すなわち、イランとモロッコ）のそれぞれについてhapFLK29（補足注5および補足図9、10、11参照）で検定し、両ケースで検出される選択下にあると考えられる共通領域に焦点を当てた。 グループのサンプルサイズ（n = 13-20）は、本方法の要件に適合していた29。 hapFLKで見つかった一貫した選択のサインが、各属の対応する世界パネルセットにも存在するかどうかを目視で確認したが、これらのグループは多品種構成であるため統計的検定に含めなかった。 最後に、層別FDR法を用いて、OvisとCapraの間で共有される選択のシグナルを探した。 その戦略は補足図4に描かれている。

各属の野生グループとイラン、モロッコの各グループを対比するためにhapFLK検定を実施した。 変種の1％のランダムサブセットを用いて、グループのペア間のレイノルド遺伝距離64から近親行列を算出した。 推定された集団樹は、neighbor-joining アルゴリズムを用いて構築された。 各SNPについて、選択的掃引の検出力を高めるためにハプロタイプ情報を組み込んだhapFLK検定を実施した。 各SNPについて、hapFLK統計量は、キンシップ行列によって推定される中立モデルに対するハプロタイプ頻度の偏差を算出した65。 連鎖不平衡の情報を利用するために、hapFLKはScheet and Stephensの多座遺伝子型に対する多点モデル66を使用しており、位相のないデータに適合させることが可能である。 このモデルの主な用途の一つは、位相推定を行うことである（fastPHASEソフトウェア66）。 我々の解析では、このモデルは位相のないデータで学習させたため、位相の不確実性を考慮した解析となっています。 この方法は、隠れマルコフモデルを用いて、染色体に沿った局所的なハプロタイプを指定されたクラスタ数K（25に設定）に再グループ化するために用いられた。

各属について、伝統的に管理されている二つの国内グループで推定的に選択されている共通領域を特定するために、先の二つのhapFLK分析を結合した。 それぞれの解析について、hapFLKスコアをχ2分布にフィットさせてp値を求めた（スクリプトはhttps://forge-dga.jouy.inra.fr/projects/hapflk/documentsで入手可能）。 野生グループと各家畜グループとの2つの対比の結果をStoufferの方法67で結合し、野生と家畜の比較のための単一のp値を得た。 最後に、FDRフレームワーク68をSNPsの全セットに適用し、結合されたp値をq値に変換した。 4804>

OvisとCapraの間で選択シグナルが共有されているかどうかを調べるために、まず2つの参照ゲノムのクロスアラインメントを用いて、相同なセグメントを同定することにした。 その後、層別FDRの枠組みを適用した69。 このアプローチは、遺伝データの事前情報を与えられた検定には固有の層別化があるという事実に基づく69。なぜなら、真の対立仮説の基礎となる分布は、様々なゲノム領域の異なる動態に従って異なり、p値の異なる分布につながる可能性があるからである。 このため、FDR調整後のp値（すなわちq値）を異なる層について個別に求める必要がある。 そこで、各属において、他属で検出された領域と相同な領域（共有層と呼ぶ）とそれ以外の領域（一般層と呼ぶ）を分離し、収束を探した。 2つの層それぞれについてp値を別々に抽出し、FDRの枠組みでq値を算出した。 これらの層別q値は、統計的有意性（<10-2）を考慮した最終量であり、選択中のSNPを検出し、対応するゲノム領域にマージした。

両属の野生動物と家畜を区別する選択の収束的なサインを検証するために、一方の属のq値（0.2から0.002まで変化させた）に適用した有意閾値と、他方の属の共有層でSNPが選択される推定確率の関係をStorey et al.70 を使って検討した。 の手法で推定した。 q値に適用する閾値の減少（ストリンジェンシーの増加）に伴って推定確率が増加するのは、ある属で有意な領域ほど、他の属でも有意なSNPが見つかる可能性が高いことを示している。

国内3グループ間で一貫性のない選択シグナルをフィルターで除去するようにした。 検出された各領域について、配列間の同一性パーセントに基づきNeighbor-Joiningツリーでクラスタリングした個体の位相付きハプロタイプを使用した。 hapFLKで検出された選択シグナルが、家畜の選択の緩和を示すのか、あるいは正の選択を示すのかを推測するために、野生グループと家畜グループの間で選択されていると考えられる各領域の塩基多様性の差(π)を推定した。 この差をΔπ指数として表し、各ゲノム領域について、野生群のπとイランおよびモロッコの家畜群のπの平均値との差から、全ゲノムについて計算した両群のπの差を差し引いた値とした：

$Delta \pi = \left( {pi _{rm wilds}}} ) – ЪЪЪЪЪЪ \γright)_{{mathrm{operatorname{genomic-region}}} – ʕ-̫͡-ʔʔ – ʕ-̫͡-ʔ-ʔ \right)_{{mathrm}}{whole-genome}}$

負の値は、2つの国内グループの平均と比較して野生グループの塩基多様性が低いことを示し、これらの最後のグループの選択の緩和、国内での多様化選択または野生での正の選択を示すとみなされるであろう。 反対に、正の値は、家畜群に生じた方向性のある正の選択、あるいは安定化する選択を示すと考えられる。 また、検出された選択的掃引がΔπ指数で与えられた示唆と一致するかどうかを、ハプロタイプ・クラスターを使って各領域で手動で検証した

機能的解釈は以下のように行った。 選択中の各領域について、機能的役割を特定するためにその領域＋両側50kb、遺伝子の上流・下流5kbを考慮し、これらの座標間の重複を評価して、目的の遺伝子を保持するようにした。 最後に、ある検出された領域と遺伝子の位置が重なっている場合に、その遺伝子と関係があると考えた。 次に、トップシグナルに最も近い遺伝子、すなわち領域内で最も低いq値の位置を考慮し、どの遺伝子が最も選択の対象になりやすいかを評価した。 遺伝子は、Uniprot（http://www.uniprot.org/）を用いて、「Biological Process」カテゴリー（すなわち、GO:0008150）の30の子項目（すなわち、その項の直下）に対する関与を考慮し、機能的注釈を行った。 哺乳類機能に関与しない3つのターム（すなわち、GO:0006791硫黄利用、GO:0006794リン利用、GO:0015976炭素利用）は考慮しなかったため、33カテゴリー中30カテゴリーについて各遺伝子に対応する全てのGOタームを検索した（補足データ4）。 GOカテゴリーに含まれる遺伝子の分布を比較するために、2つのχ2検定を行った。すなわち、（i）属特異的領域から選択された遺伝子と相同領域からの遺伝子、（ii）選択された全遺伝子とSwiss-ProtのGOタームに関連する18,689個のヒト遺伝子を比較した。 また、遺伝子の機能を家畜の文脈で解釈するために、その表現型効果に関する文献情報も入手した。

最後に、先に検出した領域内で、野生グループと家畜グループの間で最も差異が大きいSNPを見つけるために、FLK統計量を用いた。 hapFLKについては、キンシップ行列で推定した中立モデルに対する単一マーカーの対立遺伝子頻度の偏差を表している65。 2つの解析のスコアをχ2分布に当てはめ、得られたp値を組み合わせるという手順は、hapFLKの検定に用いたのと同じであった。 しかし、p値の分布が一様でないため、FDRの枠組みを適用することができず、p値<10-4を示すhapFLKで検出された領域内のSNPsを選択することにした。 これらのSNPについては、OvisについてはEnsembl v74 sheep OARv3.1ゲノムアノテーション（http://www.ensembl.org/Ovis_aries/Tools/VEP）から、CapraについてはNCBI eukaryotic genome annotation pipelineによって作られたgoat CHIR1.0 genome annotationから生成したVarious Effect Predictor（VEP）注釈71を用いた（https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/)。 SNPは遺伝子間、上流と下流（UTRを含む）、イントロンとエキソン位置に分類された。 FLK p値<10-4のSNPの分布と、選択シグネチャーの検出に用いた全SNPの分布の違いをχ2-検定で調べた。

データの利用

これらの解析に用いたOvis73サンプルとCapra72サンプルについて作成したシーケンスとメタデータデータは一般に公開されています。 一般的な情報とすべてのvcfファイルはEnsemblのウェブサイト（http://projects.ensembl.org/nextgen/）で見ることができる。 O. orientalisとC. aegagrusのすべてのFastqファイル、Bamファイル、de novoアセンブリは、European Nucleotide Archive (https://www.ebi.ac.uk/ena) でNextgenプロジェクトのアクセッションコード (PRJEB7436) で見ることができる。

コメントを残す

メールアドレスが公開されることはありません。