サンプル収集と全ゲノム配列決定
アジア全域の家畜用バクトリアラクタ105頭を対象に、全ゲノム配列決定とサンプル収集、サンプル収集、サンプル収集、全ゲノム配列決定、サンプリングを行った。 本研究では、MG州ゴビアルタイ地域から19頭の野生ラクダ、イランから4頭のヒトコブラクダを収集した（補足図）。 1および補足表1）。 家畜用ラクダは、中国内MG（IMG）、新疆（XJ）、青海から55頭、MGから28頭、KAZAから6頭、ロシア（RUS）から10頭、IRANから6頭と、できるだけ多くの主要地域をカバーするように選択された。 中国とMGでは家畜であるバクトリアンキャメルが広く利用されているため、様々な地方品種が形成されており、その中から代表的な8品種を選びました。
DNA抽出後、個々のゲノムを平均13×カバレッジで配列決定した（補足図2、補足表2）。
DNA抽出後、個々のゲノムを平均13×カバレッジで配列決定した（補足図2、補足表2）。配列リードは、バクトリアンキャメルの以前のゲノムアセンブリ29とアライメントし、バリアントコールした。 その結果、1,383万個の一塩基多型（SNP）と141万個の低インデルが同定された（補足図3）。 また、生SNPのトランスバージョン比（2.29）は、ドロメダリーで報告されている値（2.31-2.34）31よりも低かったが、フィルタリング後には2.44となり、SNPの品質向上が示唆された。 変異体の機能アノテーションでは、約63.10%が遺伝子間、33.62%がイントロン、0.94%がエキソンであった（補足表3）。 家畜のバクトリアンクタール、野生のバクトリアンクタール、ドロメダリーでそれぞれ1373万、639万、1055万個のバリアントが確認された。 ドロメダリーはバクトリアンラクダの両種と系統的により分岐しているが、現存の野生バクトリアンラクダで観察される遺伝子変異の大幅な減少やドロメダリーと家畜バクトリアンラクダ間の遺伝子流動により、家畜バクトリアンラクダの方が野生バクトリアンラクダより多くの変異体を共有している（66.73％）ことがわかった（補足図4）。 家畜ラクダでは、東アジア集団で1268万個、中央アジア集団で1161万個の変異が確認された（補足図4）。 東アジアからサンプリングした家畜ラクダは中央アジアのものよりも多かったが、各集団に私有するバリアント数は、両地域間で大きな偏りは見られなかった（P値＝0.77、両側t検定、補足表4）。
遺伝子の多様性と分化
異なる集団間でのより詳しい比較のために、まず残った他者と近い遺伝関係を示す14個体を取り除いた（補足表5）。 ドロメダリーのペアワイズヌクレオチド多様性π（図1a）は、全地域のバクトリアンラクダのそれ（0.88 × 10-3-1.11 × 10-3；付表6）より有意に高く、これまでの個別ゲノムに基づく異型接合度の推定値と対照的であった10． 重要な理由の一つは、中央アジアにおけるドロメダリーとバクトリアンラクダの交雑の慣習であると考えられる30。 バクトリアンラクダのうち、野生集団は家畜集団と比較して最も低いπ（0.88×10-3）を示した（図1aおよび補足表6）。 この結果は、イヌ24、ブタ27、ウサギ32など、通常、野生動物の方が家畜よりも遺伝的多様性が高いという多くの事例に反するが、ウマ33のように個体数が極めて少ない絶滅危惧野生動物の場合には、そのようなことが起こるのであろう。 また、中央アジアに生息する家畜集団は、東アジアに生息する家畜集団（0.95 × 10-3-1.02 × 10-3；図1a）に比べ、一般に高い多様性を示した（1.03 × 10-3-1.11 × 10-3）。 この傾向は、Wattersonのθ（補足図5）でも確認された。 ここでも中央アジアのドロメダリーとの交雑がこの地域の多様性の高さを説明していると考えられる。
a Nucleotide diversity π. boxplotはゲノム全体の10kb-スライドウィンドウ2.0×105個のπを示す． 左側には各集団の地理的起源とサンプルサイズを、右側にはπの平均値を示している。 MG, IMG, XJについては、複数の地方品種をサンプルとした。 遺伝的関係が密接な個体は削除した。 箱ひげ図の要素は次のように定義される：中心線は中央値、箱の限界は第3四分位と第1四分位、ひげは1.5×四分位範囲。 ヒートマップは2.0×105個の10kbスライディングウィンドウの平均値を示す。
次にWeirのFstによってラクダ集団間の一対の遺伝距離を測定した（図1b）． その結果、既知の系統樹とよく一致し、ドロメダリーはバクトリアンキャメルとのFstが最も高く（0.54-0.64）、野生バクトリアンキャメルは家畜とのFstが2番目に高い（0.27-0.31）ことが示された。 家畜ラクダ間の分化はかなり低く、最近の単一起源であることと一致する。 興味深いことに、国産ラクダのうち、イラン産のラクダは他のすべてのラクダと最も大きな分岐を示した（0.05-0.06）。 この集団分化を検証するために、全個体についてペアワイズIDバイステート(IBS)行列に基づく近傍結合(NJ)樹を構築した(補足図6)。
集団ゲノム推定で問題となりうるのは参照バイアスであり、単一の参照ゲノムを用いると、それと大きく異なる個体のバリアントコール効率が低くなることである34。 このバイアスを調べるために、シークエンス深さを共変数として、3種の間でバリアントの欠損数を比較した（補足表7）。 分散分析の結果、家畜のバクトリアラクダは野生のラクダと有意差はなかったが（P値=0.50）、確かにドロメダリーより欠損数が少なかった（P値=4.38 × 10-3）。 このバイアスが推定値に与える影響を評価するため、コーディング配列は種間で不変である可能性が高いため、同義語SNPのみを用いて遺伝的多様性とFstを再計算した（補足図7）。
混血を伴う集団構造
混血を伴う集団構造全体を明らかにするために、連鎖不平衡が大きく、機能的効果が期待できるSNPを除去してプルーニングを行った。 プルーニングしたサブセットに基づく多次元尺度法（MDS）分析では、フルセットと同様の結果が再現された（図2a、補足図8）。 予想通り、ドロメダリーと野生のバクテリオン・ラクダルはそれぞれ第1、第2座標で分離できた。
a Multidimensional scaling (MDS) plot with coordinate 1-4 (C1-C4). b Admixture analysis assume different number of ancestry K. (混血解析)． c 移住イベントの仮定mを変えたツリーミックス解析。移住の重みは、ドナー集団から受け取った祖先の割合である
異なる祖先割合を推定するために、Admixture35を用いて、先祖集団Kを仮定した集団構造分析を行った（図2b）。 クロスバリデーションの結果、K＝3が最適であることが支持され（補足図9）、ドロメダリー、野生バクトリアンキャメル、家畜バクトリアンキャメルに明確に分かれていることが示された。 イラン産ドロメダリーに家畜用ラクダが導入されたことは、少なくとも1頭のドロメダリーで確認された。 したがって、イラン、カザフスタン、ロシアなどの中央アジアのラクダ集団には、ドロメダリーの祖先が多く、その割合は1〜10％と推定される。 さらに、いくつかの野生個体には家畜であるバクトリアンラクダの祖先が認められ、その割合は7〜15％であった。 これは祖先の多型に起因する可能性もあるが、mtDNAs36やY染色体15で観察され、野生種の遺伝的独自性を脅かすと指摘された内生的交雑に起因する可能性もある。 意外なことに、野生ラクダと生息地が近いモンゴルの集団でさえ、家畜集団の祖先に野生ラクダはほとんど寄与していないのである。
混血を伴う集団構造を調べるもう一つの方法として、TreeMix37を用いてラクダの集団樹を推定した(図2c)。 移動経路を加えない場合、樹形はやはりイランが国産バクトリアラクダの中で最初に分離した集団であることを示していた（補足図10）。 移動トラックを増やすと（m=1-3）、モデルの適合度が大幅に向上し（補足図11）、KAZA、RUS、IRANを含む中央アジアのドロメダリーから家畜ラクダへの遺伝子フローが、移動重み4〜9％で確認された（補足表8）。 また、移動経路はIRANへのドロメダリー支流の終点を指しているが、KAZAとRUSへのドロメダリー支流の中点を指していることも特筆される（図2c）。 このことは、イラン産のヒトコブラクダに関連するゴースト集団が、KAZAやRUSの祖先に貢献したことを示唆しているのかもしれない。 さらに移動トラック（m = 4）を追加すると、モデルの適合度が引き続き向上し、ドロメダリーのXJ種への移動が示唆された（図2c）。 TreeMixは野生ラクダと家畜ラクダの間に移動の強いシグナルを検出しなかったが、残基は野生種と東アジア種の間に中程度の混血を示した(補足図11)。 次に、これらの混血事象の統計的な有意性を評価するために、あまりパラメータ化されていない3・4個体群（F3/F4）検定38を使用した。 ここでもF3検定により、KAZA、RUS、IRANにおけるドロメダリーとバクトリアンキャメルの混血が強く支持された(Supplement Table 9)。 また、より感度の高いF4検定では、中央アジアのドロメダリーとバクトリアンクタールの混血の程度が、東アジアのそれと比べて有意に高いことが確認された（補表10）。
内生化の除去による中央アジア起源の証拠
東アジアと中央アジアは考古学的証拠に基づくバクトリアンラクダの家畜化地域の2つの選択肢であったが，最も可能性の高いものは未解決であった． イランの集団は他の家畜集団と最も大きな遺伝的差異を示したが、中央アジアにドロメダリーとバクトリアンラクダの混血が存在すれば、起源推定への支持を弱めることになる。 この影響を軽減するために、「BABA/ABBA」検定39 を用いて、バクトリアンラクダのゲノムからドロメダリーの導入セグメントを除去することを試みた。 東アジアと中央アジアの集団をグループ分けし、両集団とドロメダリーの対立遺伝子の共有を比較した(図3a)。 ドロメダリー1頭にはバクトリアンラクダの祖先が、野生ラクダ3頭には家畜ラクダの祖先が含まれ（図2b）、交絡因子となるため、解析ではこの4個体を除外した。 導入されたセグメントを見つけるために、PattersonのDのロバスト版である統計量fdを用い40、ジャックナイフ法を用いてZスコア＝2という厳密な有意水準を適用した（補足図12）。 その結果、中央アジア集団では11,711個（Zスコア> 2）、東アジア集団では3891個（Zスコア< -2）であり、予想通り導入シグナルと考えられる個が多く存在していることが判明した。 そこで、残りのセグメントをもとにAdmixture解析を行ったところ、ドロメダリーの導入が効果的に抑制されていることが確認された（補足図13）。 内殖を除去した後にペアワイズFstを再計算したところ、やはりIRANが国内全集団の中で最も分化している（0.04-0.06）ことが示された（図3b）。 さらに集団の系統を知るために、ペアワイズFstに基づいてNJ木を再構築し、ブートストラップ検定を行った（図3b）。 その結果、国内の全バクトリアン集団の中で最初に分離したのはIRANであり、次いでKAZA、RUSであることが確認された。 また、系統樹に基づくベイズ型二元マルコフ連鎖モンテカルロ（MCMC）解析では、中央アジアが家畜用ラクダの祖先地域であり（確率＝99.78%）、その後の分散経路も中央アジアから東アジアであることを強く支持した（補足図14）。
a BABA/ABBA分析でドメダリーの内生化を家畜用ラクダへの導入を確認したもの。 この導入に注目するため、バクトリアンラクダの祖先を持つドロメダリー1頭と、家畜の祖先を持つ野生ラクダ3頭を取り除いた。 各樹形において、有意なfd（｜Z-score｜ > 2）を示した100kbセグメントの数を示す。 b 近接したセグメントを除去した後の集団のNJ（Neighbor-Joining）樹。 ヒートマップは5.1×104個の10kb-sliding windowの平均ペアワイズFstを表す。 NJ樹のブートストラップ値は、5,000個の10 kbの窓を100回ランダムにサンプリングして算出した。 母集団は異なる色で表し、Genbankからの配列はドットで示した。
独立した証拠として、我々は、本研究で配列決定した128サンプルと、Genbankから利用できる39の追加サンプルに基づいて、完全長mtDNAの最尤樹も再構成した（図3cと補足表11）。 mtDNAの導入は容易に特定でき、ツリーから除外することができた。 例えば、KAZAとRUSから新たに配列決定された2つのmtDNAは、ドロメダリーとクラスター化した。 家畜ラクダのクレード内では、異なる地域のラクダがほとんど混在していたが、イランからの2つのmtDNAが家畜集団の最も基部の枝を形成していた（図3c）。
バクトリアンラクダの人口動態
我々は歴史上のラクダの人口動態を推測するために、いくつかのパラメトリックモデリング解析を行った。 先行研究10と同様に、pairwise sequentially Markovian coalescent（PSMC）モデル41に基づくバクトリアラクダの長期軌跡は、100万年前より早く祖先ラクダの有効集団サイズが驚異的に減少することを明らかにした（補足図15）。
ラクダ集団間の分岐時期をより詳細に調べるため、一般化系統樹合体サンプラー（G-PhoCS）42を使用した。 G-PhoCSは、ラクダ集団の系統が与えられたとき、各集団の個々のゲノム中のリンクしていない中立遺伝子座に基づいて、突然変異でスケールした集団サイズと集団分岐時間を推定することができた。 モデルの複雑さを軽減するために、ドロメダリー、野生のバクトリアンキャメル、および家畜バクトリアンキャメルの代表的な3つの集団（IRAN、KAZA、MG）だけを対象とした（補足図16、補足表12）。 図3bによると、IRANとKAZAは中央アジアで最初に分離した集団であり、MGの分離は中央アジアから東アジアへの分散を示す可能性がある。 年代はTimetreeデータベース43により、バクトリアラクダとドロメダリーの分岐を573万年と仮定し、較正を行った。 移動帯を組み込まない場合、すべてのパラメータ推定値の収束は容易に達成できた（補足図17、補足表13）。 PSMCの結果と同様に、祖先から現代の集団まで、有効人口規模は概して減少していた（図4）。 野生ラクダと家畜ラクダの分岐時期は0.43百万年前（95％信頼区間：0.13-0.73Mya、図4）と推定され、mtDNAに基づく時期（0.714または1.1Mya9）よりやや遅かった。 国内集団のうち、イランは約445万年前（95%CI: 0.07-17.6 Kya）に他と分離し、その後、中央アジアと東アジアの集団は約240万年前（95%CI: 0.01-7.84 Kya、図4）に分離していることがわかった。
突然変異スケール有効人口サイズ θの変化は熱色により表現されています。 年単位の時間は、ドロメダリーとバクトリアンキャメルの分岐時間で較正した。 95%信頼区間は時間軸の棒グラフで示した。 赤と青の棒はそれぞれIRAN-MGとKAZA-MGの分岐を示す。
遺伝子流動を考慮し、ドロメダリーからバクトリアンキャメル集団への移動帯も導入してみた（補足図16、補足表12）。 イラン産のドロメダリーからIRANへの移動バンドとゴースト集団からKAZAへの移動バンドを導入した場合のみ、推定値が収束した（補足図18、補足表13）。 野生ラクダと家畜ラクダの分岐時間は、移動モデルによって変化しなかったが（0.46Mya, 95% CI: 0.24-0.71Mya） 、家畜集団の最初の分岐時間（0.19Mya, 95% CI: 0.08-0.31Mya ）は、既知の家畜化の歴史（11.5 Kya44）をはるかに超えており、非現実的なものとなっている。 また、移動率はイランへの移動帯で0.27％、カザへの移動帯で0.16％と、Admixtureによる推定値（1〜10％）よりはるかに低い値しか求められなかった。 G-PhoCSが想定する一定率の連続移動モデルよりも複雑な混血の歴史が、この推定値の低さの理由であると考えられる。