Shared polymorphism is abundant between A Arabidopsis and its relatives.This case is a long term balancing selects in the American American Genome Biology. 7878>
A. thalianaの80アクセッションからなる集団において、4,902,039のSNP（119,146,348サイト中）が存在し、そのうち2,044,731はマイナーアレル頻度が> 0.05 であった。 C. rubella集団では、22のC. rubellaアクセッション（Additional file 1: Table S1、21の公開アクセッションと本研究で配列決定した1アクセッションを含む）から、C. rubella参照ゲノムに対してSNPをコールし、2,149,643SNP（134,834,574部位から）を特定した、このうちMAF > 0.05 のものは 1,240,547 件であった。 2種間で共有される多型（特定のオーソログ部位で同一の対立遺伝子対として定義）を同定するために、まず2種間のオーソログ遺伝子対のセットを構築した。 オーソログ遺伝子が保存されていることを保証するために、A. thaliana と C. rubella の参照ゲノムに加え、A. thaliana の近縁種である Arabidopsis lyrata を含めた。 その結果、16,047組のオルソログ遺伝子が得られたが、3つのリファレンスのいずれかにタンデム重複がある33組を削除し、最終的にA. thalianaとC. rubellaのオルソログ遺伝子16,014組を得て、さらに解析した。
A. thalianaのオルソログ遺伝子16,014個のゲノム領域は39,275,210 bp、同様にC. rubellaのそれは40,936,262 bpにわたっている。 これらの領域には3,889,495の固定的な違いがあり、この高い比率（約10％）は、2つの種の長い分岐時間（〜8 MYA）と一致する。 これらの領域では、A. thalianaでは1,122,845のバイアレリックサイト（MAF > 0.05の426,123）、C. rubellaでは452,116バイアレリックサイト（MAF > 0.05の279,780）を発見しました。 これらの多型部位のうち、19,732のオーソログ部位が両種において多型であり、そのうち8535は同じ対立遺伝子ペア（共有SNP ）を共有していた（追加ファイル1：表S2）。
非コード領域配列と比較して、コード領域配列は保存度が高く、分岐が激しい2種間で強固なアライメントが得られたため、まずコード領域のshSNPsに注目した。 両種ともMAF > 0.05はSNPの信頼性を保証し，長期的な平衡選択下にある部位では中間頻度のアレルの過剰が予想されるため，これを考慮する必要があった。 1007遺伝子のコーディング領域に1503個のshSNPを発見した。
ジェノタイピングとマッピングエラーを避けるために、1503個のshSNPにさらにフィルタリングを適用した。 A. thalianaのSNP行列をダウンロードしたため、フィルタリングはC. rubellaのSNPデータのみに適用された。 ゲノム重複による偽SNPを避けるため、C. rubellaの50bp領域すべてのマッピング可能性を評価し、一意にマッピング可能な領域内にあるサイトのみをその後の解析に残した。 その結果、580サイトが残った。 最後に、SNPコーリングツールでマークされた低品質の部位を除去した後、433の遺伝子に信頼性の高い共有コーディングSNPを546個得た。 フィルタリングの詳細については、「方法」のセクションに記載されており、その様子を図2に示す（
Pipeline of the SNP filtering process to identify candidate TSP sites
2種の人口動態履歴
豊富な共有多型から本物のTSP信号を検出するには、2種の人口動態履歴を完全に理解することが必要である。 ジョイントサイト周波数スペクトル(joint SFS)は、多様な生物の人口動態史を研究するために広く用いられている 。 そこで、まず、A. thalianaとC. rubellaの基準ゲノムのアラインメントから、16,014個のオーソログについて、4倍退化部位を抽出した。 そして、遺伝子流動を伴わない基本モデル（M1、図3）と、両属間の古代遺伝子流動を取り入れたモデル（M2、図3）で、fastsimcoal2を用いて合体シミュレーションを行ったところ、2,011,573サイトが得られた（詳細は「方法」参照）。 異なる属に属し、染色体数が異なる種（5本と8本）は、最近の導入の可能性が極めて低いため、2種間の古代の遺伝子フローのみを考慮した。 また、両属とも染色体数が8本ではなく5本の種はA. thalianaのみであることから、A. thalianaが他のArabidopsis属から分離する前の古代の遺伝子フローに限定した。 各モデルにおいて、2つの属の分岐時間を8MYA（800万世代前）とし、自然突然変異率を7×10-9/bp/世代と仮定した。 A. thalianaは6MYA頃にシロイヌナズナ属から分岐した後に個体数が減少し、C. rubellaはC. grandifloraからの種分化に伴ってごく最近ボトルネックを経験した。 我々は、fastsimcoal2に実装されている複合尤度法を用いた合体シミュレーションを行い、抽出された2,011,573個の種を越えた4倍体縮退部位から計算された2種のジョイントSFSに両モデルを適合させた。 Excoffier et al.と同様に、Akaikeの情報量基準（AIC）とAkaikeの証拠能力（w）を用いて、2つのモデルを比較した。 その結果，古代の遺伝子フローを考慮しないモデル（M1）の方が若干適合度が高く（Max EstLhood: -682010 vs -682028），AICも低く，重みも大きかった（図3，追加ファイル2：表S3）． また、2つの尤度が近いことから、祖先の遺伝子流の影響は長い時間スケールで一掃されているはずで、モデルの品質にはほとんど寄与していないことがわかる
Demographic parameter estimates for two models of the divergence of the two species
Model M1 では、現在の A. N e. thalianaの現在のN eは519,000、95％信頼区間（CI）は486,368-527,574で、大きな祖先集団（〜2,230,000、95％CI = 1,085,330-4,876,051）から、〜 5.84 MYA (95% CI = 5.27-6.70) で残りのArabidopsis genusと分離していることが明らかになった。 C. rubellaは〜0.40MYA (95% CI = 321,998-500,317) に、N eが〜4,037,000 (95% CI = 2,076,868-5,165,614) で現在のN eが〜129,000 (95% CI = 126,383-157,779) と大きな祖先集団から進化を遂げた。 この2つの属はN e = ~ 4,930,000 (95% CI = 4,560,931-4,969,696) の祖先集団から分岐したのである。 遺伝子フローを考慮したモデルM2では、シロイヌナズナ属の祖先のN eが大きく（〜 3,270,000, 95% CI = 797,016-4,342,346 ）、Capsella属のN eが小さく（〜 1,972,000, 95% CI = 2,126,346-6,248,003 ）なった以外は同様のパラメータ推定値が得られた。 また、Capsella属からArabidopsis属への遺伝子フローは、逆方向よりも強いと推定された（世代あたりの移動率；1 × 10-8, 95% CI = 4.0 × 10-15-1.1 × 10-6 vs 7 × 10-14, 95% CI = 5.7 × 10-15-6.1 × 10-5）、ただしどちらも弱い（詳細は、追加ファイル2：表S3参照）。
2種間の種間多型は平衡選択下になければならない
種間多型は中立である可能性があり、その確率は特定の人口動態パラメータを与えて近似することが可能である。 ヒトとチンパンジーのTSPの研究と同様に、中立進化の下では、我々の系では以下の場合にのみ共有多型が子孫によって同一となる。 (1) 少なくとも2つのA. thaliana系統と2つのC. rubella系統がA. thalianaとC. rubellaの分裂前に合体していない、 (2) 同じ対立遺伝子を持つ系統が異なる対立遺伝子を持つ系統より先に合体している。 この確率は主に条件(1)によって決まり、合体理論に基づいて次のように近似できる：
ここでTは2属間の分岐時間、N A/N CはA. falumの集団サイズであり、P={e}^{-sec]であった。 thaliana/C.rubellaの集団サイズをそれぞれ示す。 モデルM1での推定によれば、個体数の変化を考慮すると、この子孫による同一性の確率は10-9のオーダーになる。 遺伝的領域で2つの種の間に< 39,275,210の配列部位があることから、遺伝的ドリフトのみによって中立的TSPの総数が< 1になると考えられる。
我々のモデルではランダム交配を仮定したが、両種は自家発生であり種内の集団構造はおそらく存在すると考えられる。 それでも、ゲノムの同じ領域で偶然に両種に深い合体が起こることを必要とするので、最近の人口統計学的イベントは比較的影響が少ないはずである。 また、以前の研究で示されたように、現代人における深い集団構造も確率にほとんど影響を与えないはずである。 本研究では、両種とも外交配が優勢であった歴史を持つ。 A. thalianaが外交配から自殖に移行したのはわずか100万年前であり、C. rubellaが移行したのはもっと最近である 。 また、自殖の種であっても、局所的な個体群の交雑率は14.5%という高さである。
釣り合い選択下の種を超えた多型の同定
TSPは、長期の釣り合い選択下にある領域が種ごとではなく、対立遺伝子ごとに集まるため、中立変異と区別することができる 。 そこで、コーディング領域で信頼できる共有SNPを持つ433の候補遺伝子に注目し、遺伝子領域でMAF > 0.05の共有バイ・アレリックSNPを覆うハプロタイプを調べた。
TSPの信号を持つ各セグメントの長さを推定するために、以前に導き出した、組み換え率に大きく依存する公式を用いた。 合体という観点から見ると、このようなセグメントは、同じ対立遺伝子クラスからの系統がすべて祖先集団の最も新しい共通祖先に合体するまでは、組換えによって分割されることはない 。 両種について3.6 cM/Mbの組換え率を採用すると、セグメントの長さは極めて短く、つまり理論的には数塩基対に過ぎないことがわかった。 両種が最近それぞれの外来種から発生し、過去の有効組換え率がはるかに高かった可能性を考えると、予想される長さはさらに短くなる可能性がある。 この推定値は、我々のシステムにおける中立的な状況下では、組換えによる中断のないセグメントを発見することは非常に困難であることを示唆している。 しかし、バランス選択が存在する場合、選択によって周辺領域の組換えを抑制することができる 。 そのため、セグメントの長さは、中立的なモデルで推定される理論値よりも長くなるはずである。 そこで、100bpのウィンドウサイズと1bpのステップサイズを用いて遺伝子領域をスキャンした。
433の候補遺伝子において、975の共有バイアレルSNPs（エクソンおよびイントロンを含む、MAF > 0.05）を検出した。 先行研究と同様に、次に、適格なウィンドウ（長さの95%以上で整列、詳細は「方法」参照）の中から、両種において強い連鎖不平衡（r 2 > 0.5）である975個のSNPのうち、少なくとも2つを覆うウィンドウを探して、対立遺伝子ツリーを特定した。 これらの制限により、偽陽性を大幅に減少させ、対立遺伝子が存在する場合には、高い解像度で対立遺伝子の木を得ることができる。 最後に、AT1G35220、AT2G16570、AT4G29360、AT5G38460、AT5G44000の5つの遺伝子から、10部位を含むウィンドウを、長期バランス選択下のTSP候補として特定した（追加ファイル3：図S1）。 今回見つかった5つのオーソログ遺伝子はいずれもコピー数変動（CNV）と相関がなく、2種のリファレンスとそれぞれ比較したところ、いずれも1件しかヒットしなかった（詳細は「方法」参照）。
同定した領域を検証するために、まず各集団から同定した領域のすべてのハプロタイプを求め、ハプロタイプごとに代表アクセッションを再シークエンスした（プライマーは「付加ファイル1：表S4」参照）。 予想通り、5つの遺伝子におけるTSP候補部位はすべて検証され、候補領域における2種の配列は、種ではなくアレルによってクラスタリングされた(図4)。 遺伝子AT1G35220では、2つのTSP候補部位はイントロン領域で完全な連鎖不平衡にあった。この領域はバランス選択の標的になっているか、検出されていないコーディングTSP部位にリンクしているのかもしれない。
5つの遺伝子のすべての候補領域は、種樹ではなく、対立遺伝子樹を作り出す
それぞれの領域のハプロタイプは対立遺伝子によってクラスター化したが、その結果、対立遺伝子によってクラスター化した。 AT2G16570を除いて、2つの種の間でハプロタイプの共有が検出されることはほとんどなかった（Col-0はいくつかのC. rubellaのアクセッションとハプロタイプを共有していた。） これは、分岐時間が長いことを考慮すれば驚くべきことではなく、広範なハプロタイプ共有は通常もっと短い時間スケールで現れ、近縁種間の最近の導入などのイベントによって誘発される。
Neutral simulation studies validate the five candidate genes
観察した窓が中立進化下でランダムに生成されて偽陽性となるかどうかを確かめるため、fastsimcoal2（追加ファイル4：テキストS1）を用いて推定人口動態パラメータに基づいて追加のシミュレーションを行った。 中立的な突然変異とは別に、遺伝子流は共有SNPsをもたらす可能性がある。 そこで、モデルM1（遺伝子フローなし）とM2（古代の遺伝子フローあり）の両方でシミュレーションを行ったが、人口統計学的解析の結果、M1の方がわずかにデータに適合していた。 両シミュレーションとも、変異のクラスごとに変異率の不均一性を考慮し、特にCpG部位での変異率が高く、偽陽性をもたらす可能性があることを考慮した（追加ファイル1：表S5、追加ファイル4：テキストS1）。 Fastsimcoal2を用いて、それぞれのモデルで100bpの中性セグメントを100万回生成し、TSPを探索するように、2つ以上の共有SNPを持つものを探し、アリルでクラスタリングした。 また，中立的な共有SNPが存在するにもかかわらず，どのウィンドウも対立遺伝子ツリーを形成しなかった．なぜなら，共有SNPを持つすべてのウィンドウは，2つの種の間でより多くの固定的な差異を伴い，多様性よりも高い分岐レベルを示唆していたからである． この結果は、これらのシミュレーションされた中立的な共有SNPはTSPではなく、再発突然変異であることを示唆し、さらに重要なことは、我々が見つけた5つの遺伝子は中立進化と一致せず、それによって平衡選択下の本物のTSPであることが証明された。 最終的なTSPの部位と遺伝子を表1に示す。 さらに、前述の人口統計学的研究と合わせて、この結果は、たとえ古代の遺伝子流動があったとしても、中立進化下では、この系ではTSPがドリフトによって失われることを示唆している。
バランス選択下の遺伝子の性質
次に各種の5遺伝子の全TSP領域について塩基多様性（π）計算し、M1下の中立配列でシミュレーションしてバックグラウンド多様度レベルを決定した。 5つの遺伝子のすべての領域は、C. rubellaとA. thalianaの両方で、A. thalianaのAT5G38460を除いてバックグラウンドレベルより有意に高いπ値を示した(Wilcoxon-Mann-Whitney検定、FDR補正P < 0.05, Table 2, Additional File 3: Figure S2A)。 また、これらの遺伝子の対立遺伝子は中間頻度の傾向を示した（Wilcoxon-Mann-Whitney test, P = 0.0752/0.03474 for A. thaliana/C. rubella; Additional file 3: Figure S2B）。 しかし、バランス多型にリンクした部位の対立遺伝子頻度分布は、どの対立遺伝子頻度でもあり得る頻度均衡へのシフトを示すと予想されるため、中間頻度はバランス選択を示すが、決定的な証拠とはならない 。
本研究で長期バランス選択を受けている5遺伝子のうちの一つAT1G35220は機能不明であるがエチレン処理でタンパク質リン酸化が見られる． 中でもAT2G16570はプリンヌクレオチド生合成経路の重要な酵素で、細胞分裂、葉緑体生合成、種子発芽に重要である；AT4G29360はO-グリコシルヒドロラーゼファミリー17タンパク質で防御反応に関与している． AT5G38460はグリコシルトランスフェラーゼで、ある化合物（ドナー）から別の化合物（アクセプター）へのグリコシル基の転移を触媒し、生物的ストレスを含む多様な機能に関与する ; AT5G44000はグルタチオンS転移酵素で、通常は生物的および生物的ストレスに対する反応に関与している . どうやら、これらの遺伝子は、生物学的または生物学的ストレスへの応答（AT4G29360、AT5G38460、AT5G44000）または基本的な生化学的機能（AT2G16570）に関与する可能性がある。
予想通り、バランス選択下の遺伝子は機能的に重要で、5遺伝子のすべての相同物は緑色植物の最も新しい共通祖先にすでに存在した。 表S7（Additional file 1: Table S7）に示すように、Physcomitrella patensまでさかのぼることができるAT4G29360を除く5つの遺伝子は、緑色植物の最も基層の種、Chlamydomonas reinhardtiiでも相同性（オルソログまたはパラログ）が見いだせる。 これは予想されたことで、これらの遺伝子座は短時間のリードで特定するにはあまりに可変的であるためである。 また、S-locusはシロイヌナズナゲノムの最新のアノテーションでは存在せず、C. rubellaでは、外交配から自家交配に移行し、自己不和合性が崩壊して以来、S-locusハプロタイプは1つしか維持されていない 。 さらに、両種が自殖するようになったため、S-locusはもはや平衡淘汰の下にない。 一方、今回同定した遺伝子は、古いものではあるが、これまで包括的な研究がなされておらず、どのような遺伝子が平衡選択を受けているのかを知る手がかりになると考えられる。
Balancing selection contributed to adaptation to divergent habitats
長期の平衡選択を受けている対立遺伝子変異が生態系の多様性と関連しているかどうか、48種の生態的要因に関する分岐を調べた（追記5：表S8A）。 GPS情報がないことと、C. rubellaのサンプルサイズが小さいことから、この解析はA. thalianaのサンプルに対してのみ可能であった。 個体群構造は通常、生態系の多様化と高い相関があるため、我々の結果を混乱させる可能性がある。 このような構造は、A. thalianaとC. rubellaの種樹を観察する確率には影響しないが、我々はまず、A. thalianaのサンプルにおいて、どのTSPサイトも集団構造と相関があるかどうかを確認した。 ADMIXTURE を用いると、80 個の A. thaliana サンプルは 2 群に分類でき（Additional file 3: Figure S3; Additional file 6: Table S9）、遺伝子 AT5G38460 からの 2 部位の対立遺伝子分類のみが集団構造と有意に相関することがわかった（カイ二乗検定、FDR 補正 P < 0.05,; Additional file 1: Table S10）。 そこで、AT5G38460をその後の生態解析から除外した。
生態的分岐を完全に理解するために、最近公開された1135のA. thalianaゲノムを使用した。 まず、すべてのサンプルがその自然な生息環境を高度に代表していることを保証するために「間引き」プロセスを適用し、584サンプルを残した（「方法」参照）。 次に、各遺伝子について、2つのTSP部位に対する位相差ハプロタイプに基づいて、A. thalianaの584アクセッションを2つのグループに分類した（追加ファイル5：表S8B、C、いくつかのサンプルは位相差できなかったので削除された）。 次に、4つの遺伝子それぞれについて、48の生態的要因に関して、2つのアクセッションのグループ間の分岐を評価した。 興味深いことに、これら4つの遺伝子はすべて、ある特定の生態学的パラメータの分岐と関連していた。 特にAT1G35220とAT4G29360は、温度関連の生態的要因のほとんどに関して有意な分岐を示した（追加ファイル5：表S8 A, Wilcoxon-Mann-Whitney test, FDR-corrected P < 0.05）
次に、4遺伝子すべてについて生態的ニッチをモデル化した。 その結果、ニッチの類似性を示すWarrenのI統計で示される各遺伝子の2つのサンプル群は、100個のランダムな並べ換えよりも有意に低いニッチ同一性を観測した（1サンプルt検定、FDR補正P < 0.01； 図5a、追加ファイル5：表S8 D）。 つまり、2つの対立遺伝子群のサンプルは、有意なニッチ分岐を示す。 さらに、各遺伝子の各アレリックタイプのサンプルは、小さな局所領域に孤立するのではなく、散在していた（Additional file 3: Figure S4）。 これらの結果は、これらの遺伝子座のすべてが適応と相関していることを示唆している。
生態と発現の分岐。 a 4つの遺伝子のそれぞれについて、観測されたIスコア（I O）とシミュレーションされたIスコア（I S）で示される2種類のサンプル間の有意な生態学的分岐。 b 遺伝子AT5G44000の発現分岐。 c 左：AT5G44000について2種類のサンプルの高確率（≧0.5）でのニッチのモデル化。 右側。 異なる並べ換え戦略の下での有意性の結果（確率≧0.5のニッチについて；I O = 0.673, 100回の並べ換え）
また、A. A. の84の公表された葉組織抽出トランスクリプトームを選び、2つのTSP部位の段階的ハプロタイプに基づいて2つの対応グループ間での4遺伝子の発現差異について検討した。 Thalianaの84の葉組織抽出トランスクリプトーム（各アクシオンにつき1サンプルずつ配列決定し、発現量はマッピングされた100万フラグメントあたりのエクソン断片として測定）を選び、前回と同様の検討を行った。 そこで、AT5G44000のニッチモデリングを行い（図5c）、2つのサンプル群（503対75）の多様化を検討した。 まず、高確率（≧0.5）のニッチに限定して、AT5G44000の2つのハプロタイプ群間のニッチ同一性を比較したところ、同様の結果が得られた（図5c、追加ファイル5：表S8 D）。 サンプルサイズの不均衡が結果に影響するかどうかを確認するために、各反復において両方のセットで同じサンプルサイズ（75）に分析を制限することによって、別の並べ替え戦略を使用しました（確率> 0.5で）。 Fig.5cに示したように、実サンプル群に対して並べ替えを行った場合（シミュレーション1）、観測されたI値（0.673）は有意差を示さず（1サンプルt検定、P = 0.166）、サンプルサイズの違いにかかわらず、観測値は信頼できることが示されました。 2つの実グループを混合し、実サイズのランダムな2グループを選択した場合（シミュレーション2）、または、同じサイズ（75）のランダムな2グループを選択した場合（シミュレーション3）、観測値と順列の差は再び有意となった（1サンプルt検定、P = 1.9×10-75 for simulation 2、P = 2.6 × 10-75 for simulation 3）。 これらの結果は、AT5G44000の機能的に分化した2つのハプロタイプ群は、分岐した生態学的生息地に適応したことを示唆している
。