Improved Brassica rapa reference genome by single-molecule sequencing and chromosome conformation capture technologies
Genome assembly
ゲノムアセンブリの指針とするために、B. rapaゲノムサイズをフローサイト法により、リファレンスとして稲を用いて推定しました。 その結果、B. rapaのゲノムサイズは455 Mbと推定された（Supplementary Table S1）。 さらに、BioNanoのデータをもとに作成したコンセンサスマップの全長を計算した結果、ゲノムサイズは442.9 Mbとなった（Supplement Table S2）。 また、BioNanoのデータから作成したコンセンサスマップの全長を計算したところ、ゲノムサイズは442.9Mbであった（補足表S2）。88 Gb）、BioNanoデータの456倍（〜207.70 Gb）、Hi-Cリードの164倍（〜74.64 Gb）を用いてB. rapaゲノムを組み立てた結果、1476コンティグ、コンティグN50 1.45 Mb、全長 351.06 Mbとなりました（Table 1）。 その後、Hi-Cリードを用いて22のコンティグに不一致を検出した（Supplementary Table S4）。 これらのコンティグを削除するのではなく、矛盾する領域で分割した。例としてContig01464のデータを示す（補足図S1）。 rapa genome assemblies
BioNanoマップとmate-pair reads（BRADより、http://brassicadb.org）を用いてスキャフォールドしギャップサイズを推定した結果、スキャフォールドN50 4.44 Mbの1301スキャフォールドを得た（Table 1）。 得られたスキャフォールドを染色体上に配置するために、Hi-Cデータと改良された遺伝地図を用いてスキャフォールドのアンカーリングを行った（方法論参照）。 その結果、10本の染色体上の298.19Mbの配列が固定され、その中にはHi-Cデータによってクラスタリングされた200のスキャフォールドと遺伝地図によって割り当てられた8のスキャフォールドが含まれていた。 最終的なアセンブリはB. rapa genome v3.0と呼ばれ、合計353.14 Mbの配列で、396のギャップ（2.08 Mb）があった（Table 1）。 その結果、B. rapa genome v3.0はv1.5より長く、v2.5より短いことがわかった。 その結果、CEGMA16を用いたCore eukaryotic genes (CEGs)の検索により、アセンブリの完全性を検証しました。 その結果、248個のCEGのうち247個が完全で、1個が部分的であり、我々のアセンブリですべてのCEGを検出できることが示された（補足表S6）。 次に、B. rapaのexpressed sequence tags (EST)の配列（NCBIのdbESTからダウンロード）をマッチングしてゲノムの品質をテストしたところ、99.34%のESTが新たに組み立てられたB. rapa genome v3.0で検出されることが分かりました。
Contiguity improvement
B. rapa genome v3.0 ではギャップとコンティグサイズに関して改善されています。 B. rapaゲノムv1.5はIllumina配列から作成されましたが、アセンブリv2.5ではより多くのIlluminaリードと比較的少量のPacBio配列データが使用されました。 これら2つのアセンブリには、断片化や低い連続性による限界がある（表1）。 1分子シーケンス、オプティカルマッピング、Hi-Cテクノロジーを組み合わせることで、B. rapaゲノムv3.0は、以前の2つのアセンブリに比べて、約27倍（contig N50: 1446 Kb vs. 53 Kb, v2.5）、約31倍（contig N50: 1446 Kb vs. 46 Kb, v1.5）に連続性が改善しました（表1）また、それぞれのアセンブリのギャップサイズと量も評価されました。 v3.0では、既知のギャップ（BioNanoから122個、mate-pairデータから74個）と未知のサイズ（Hi-C scaffold接合から190個、遺伝地図接合から10個）を含めて、わずか396個のギャップがありました。 以前のアセンブリと比較して、v3.0では1Mbあたりのギャップサイズが約10倍（5.89Kb vs. 60.59Kb, v2.5）、7倍（5.89Kb vs. 40.09Kb, v1.5）改善されました（表1）。 また、v3.0はv2.5、v1.5と比較して、1Mbあたりのギャップ数が23倍（1.15 vs. 25.98, v2.5）、35倍（1.15 vs. 40.09, v1.5）少なく、それぞれ優れています（Table 1）
3 バージョンの B. B. の隣接性とスキャフォールド順序の正確さを評価するために、v3.0は、1Mbあたりのギャップの数が1倍（2.5 vs. 40.09, v1.5 ）少なく、それぞれ優れています（Table 2）。 8007＞3つのバージョンのB. rapa参照ゲノムの隣接性と足場順序の正確さを評価するために、まず、2つのヘディング白菜系統の交配から得られた2倍ハプロイド（DH）集団の同じセットのリシークエンスデータを用いて、3つのアセンブリに基づいて遺伝地図を再構築した17。 次に、遺伝地図上のビンマーカーの位置を、対応する物理地図と統合することによって評価した。 その結果、892個のbinmarkerのうち、877個（98.3%）のbinmarkerが遺伝地図にマップされた。 801個（91.3%）のbinmarkerが遺伝地図と一致し、v3.0の品質の高さが示されました（図1；Supplement Table S7）。 しかし、A05、A08、A09染色体上の76個（8.7%）のbinmarkerが遺伝地図の曖昧な位置にマップされていることに気づきました。 これらの領域は、以下の解析で説明するように、特に染色体中心領域に繰り返し配列を含んでいた。 しかし、これらの矛盾する領域はPacBioリードおよび/またはBioNanoマップによってカバーされていた。v3.0のA08染色体のデータを例として示す（補足図S2）
X軸にB. rapa genome v3.0の遺伝地図のマーカーを示し、B. rapaの物理地図のマーカーはX軸に示している。 v2.5の遺伝地図には1092個のbinmarkerが、v1.5の遺伝地図には866個のbinmarkerが存在した。 しかし、v2.5とv1.5の遺伝地図にはそれぞれ88.7%（1092個中969個）、92.3%（866個中799個）しかマッピングできなかった（補足表S7）。 その結果、v2.5のbinmarkerの15.1%（969個中166個）は、同一染色体内で遺伝的距離と物理的距離に乱れがあるbinmarker（intrachromosome）146個、異なる染色体上で遺伝的距離と物理的距離に不整合があるbinmarker（interchromosome）20個を含み、不整合なbinmarkerでした（補図S3；補表S7）。 v1.5では、10.0%（799本中80本）のbinmarkerが不一致で、そのうち71本はintrachromosomeのbinmarker、9本はinterchromosomeのbinmarkerでした（補足図S4; 補足表S7). しかし、v3.0は染色体内ビンマーカーが8.7%（877個中76個）と最も少なく、染色体間ビンマーカーは不一致であった（補足表S7）。このことから、B. rapaゲノムv3.0は以前の2つのアセンブリよりも高い連続性を有することが示された。 このように、B. rapa genome v3.0は3つのB. rapa assemblyの中で最も連続性が高く、足場の並びも良いことが示唆された。
Comparison of genome annotation
遺伝子モデルの予測および注釈は既出6と同様であった。 v3.0では合計45,985のタンパク質コード遺伝子モデルを同定し、これはゲノムアセンブリの14.74%に相当した(Table 1)。 このアセンブリでは、98.75%（45,985個中45,411個）の遺伝子が染色体上に位置し、1.25%（574個中）だけがスキャフォールド上に位置していた。 v3.0のde novoアノテーションされた遺伝子は、Brassica reference genomesの遺伝子モデル命名法の標準（http://www.brassica.info/info/genome_annotation.php）に準じて命名された。 遺伝子モデル数はv1.5（41,020個）より多く、v2.5（48,826個）より少なかった（表1）。 さらにアノテーションの質を評価するために、1440の保存植物遺伝子をベンチマークとしたBUSCO18を用いて、過去のアセンブリのアノテーションとの比較が行われた。 これらの保存された植物遺伝子の約97.7%が同定され、1.7%がv3.0に提示されたフラグメントとして検出された（Supplementary TableS11）。
シンテニック遺伝子ペアとタンデム遺伝子アレイを特定するためにSynOrths19を用いて3つのアセンブリ間でゲノムシンテニー分析を行った。 v3.0では合計2077個のタンデムアレイ（4963個のタンデム遺伝子に対応）が同定された。 v1.5でも同数のタンデムアレイ（2077アレイ、5004遺伝子に対応）が検出された。 ゲノム全体のシンテニーを評価した結果、v3.0の1539個のタンデムアレイ（3757遺伝子に対応）は、v1.5の1494個のタンデムアレイ（3670遺伝子に対応）とシンテニーであることが示された。 しかし、v2.5ではより多くのタンデムアレイ（3535アレイ、8002遺伝子）が同定された（表1）。 また、3.0、v1.5では検出されなかったsuperfluous tandem遺伝子の領域に、v2.5ではギャップが検出された（図2a）。 このギャップは、v2.5ではPacBioリードを用いたギャップクロージングによるアセンブリエラーが発生し、その結果、タンデム遺伝子のアノテーションが不正になったものと思われる。 また、ギャップがない他のタンデム遺伝子については、v3.0やv1.5の単一遺伝子が、v2.5では2つ以上の遺伝子としてアノテーションされていることが確認された（図2b）。
a BraA01000818とBraA01000819の遺伝子間に25bpのギャップ（赤い矢印で示した細い黄色のバー）があり、v2.5の無効なアノテーションであることがわかる例。b v2.5のBraA02003894とBraA02003895の遺伝子は、v3.0（BraA02g039730.3C）とv1.5（Bra020703）では単一遺伝子としてアノテーションされている。 図はGEvo (https://genomevolution.org/coge/GEvo.pl)
各タンデムアレイを単一遺伝子座とした場合、v3.0では43,099遺伝子、v2.5では44,359遺伝子、v1.5では38,093遺伝子残っていました(表1)。 次に遺伝子シンテニー解析を行ったところ、v3.0の39,858遺伝子（92.48%）がv2.5の40,442遺伝子（91.17%）、v1.5の35,464遺伝子（93.10%）と対になることが判明しました。 アノテーションされた遺伝子を初期バージョンのものと比較した結果、v3.0ではv2.5およびv1.5と比較して3241のバージョン固有の遺伝子が同定されました。 このうち、2380の遺伝子はB. rapaのmRNAリード（BRADより、http://brassicadb.org/）と一致する証拠から、2295の遺伝子は他のアブラナ科のタンパク質配列から支持された（補足表S12）。 また、v3.0に含まれるバージョン特異的遺伝子の89.10%（3214の2888）がB. rapaのmRNAデータまたは他のアブラナ科のタンパク質配列によってサポートされ、10.90%（3214の326）だけがサポートされなかった。
A new LTR-RT expansion event identified in the updated assembly
Virtexの注釈は以前報告した20と同じ方法によって行われた。 v3.0では1244のファミリーから合計235,683のTEが同定され、v2.5やv1.5と比較して562のユニークなTEファミリーが発見された。 v3.0では、TEはアセンブルされたゲノムの37.51%（134 Mb）を占め、これまでのアセンブル（32.30%、126 Mb, v2.5; 25.44%, 72 Mb, v1.5）より高い値でした2,6. 今回の新規アセンブリでは、最も豊富なTEはLTR-RTで、全長57.64 Mb、アセンブリゲノムの16.32%を占めた。 非LTR-RTリピート（LINEとSINE）は3.10%であった（補足図S5）。 26.35Mbに相当するDNAトランスポゾンを検出し、組み立てられたゲノムアセンブリの7.46%を占めた（補足図S5）。 v3.0で同定されたTEとリピートの完全なリストは、補足表S13にあります。 さらに、B. rapaゲノムv3.0では、1231個のmiRNA、1281個のtRNA、2865個のrRNA、3737個のsnRNAを同定した（補足表S19）。
今回の組み立てでは、v2.5（44Mb）およびv1.5（18Mb）と比較してより多くのLTR-RT（57Mb）の注釈付けが行われた。 v3.0では51,062個の非インタクトなLTR-RTが同定された。 さらに解析したところ、non-intact LTRの65.27% (51,602 個中 33,672 個) が10本の染色体上に存在し、non-intact LTR-RT の34.73% (51,602 個中 17,922 個) がアンカー足場上に存在することが判明しました。 同じ方法6で、v3.0では合計13,318本の無傷のLTR-RTがアノテーションされました。 しかし、v2.5では4129個、v1.5では801個の無傷のLTR-RTが見つかっただけでした6。 さらに解析の結果、インタクトなLTR-RTの18.19%（13,318本中2423本）しか10本の染色体上に存在しないのに対し、v3.0ではほとんどの（81.81%, 13,318 本中 10,895 本）インタクトなLTR-RTがアンカー足場上にあることが判明しました。0.インタクトなLTR-RTの挿入時期を既報4と同様に計算したところ、B. rapaゲノムはB. oleraceaから分岐してから3回のLTR-RT拡張の波があったことがわかった(Fig. 3)。 これらの無傷のLTR-RTの平均挿入年代は188万年前（MYA）で、挿入年代の中央値は1.59MYAであった。 さらに、v3.0ではv2.5やv1.5と比較して、異なる長さの無傷のLTR-RTが多く見つかった（補足図S6）
これらの無傷のLTR-RTにより、B. rapaゲノムにおける新しいLTR-RT拡大イベントが特定された。 0 MYA から 0.4 MYA までの 3155 個の無傷の LTR-RT 挿入イベントを、平均長 8135 bp、平均挿入日 0.2 MYA の「若い拡張」、1.0 MYA から 1.4 MYA までの 2283 個の無傷の LTR-RT 挿入イベントを、平均長 11902 bp、平均挿入日 1.2 MYAの「中位の拡張」、および 3.4 MYA の 1444 個の無傷の LTR-RT 挿入イベントとして指定した。若い拡張と古代の拡張は、これまでに同定された拡張イベントと密接に対応している。中位の拡張はB. rapaゲノムで最初に同定され、B. oleraceaの無傷のLTR-RT拡張イベントと挿入時期が類似している。 さらに、v3.0では1778個のTy1/Copia様LTR-RTと4179個のTy3/Gypsy様LTR-RTが同定され、これは以前のアセンブリで同定されたものよりもはるかに多い（v2.5では353個のTy1/Copiaと632個のTy3/Gypsy、v1.5では260個のTy1/Copiaと162個のTy3/Gypsy）（補足表 S20；補足図 S7、S8）。 一般に、Ty3/Gypsy様LTR-RTはTy1/Copia様LTR-RTよりも多かった（Supplementary Table S20）。 v2.5やv1.5と比較すると、v3.0のTy3/Gypsy様LTR-RTは5MYAから明らかに増加し（補足図S7）、Ty1/Copia様LTR-RTは2.2MYAから増加していることがわかった（補足図S8）。 系統樹から、各LTR-RT群はv3.0ではv2.5やv1.5よりもコピー数が多いことがわかった（補足表S21、S22; 補足図S9、S10）。
異なる時期（100万年前、MYA）に生まれた無傷のLTR-RTの数は、the B. の3つの集合体においてである。 ラパゲノムとB. oleraceaゲノムの3つのアセンブリにおける、異なる時期（100万年前）に生まれたLTR-RTの数。
B. rapaゲノムにおけるゲノムブロックとセントロメア
最新のアセンブリv3.0を用いてゲノムブロック間の関連性を検討しました。 B. rapa genome v3.0のゲノムブロックとセントロメアを定義するために、まずv3.0とA. thalianaのシンテニック関係に基づいて3つのサブゲノム（LF、MF1、MF2）を構築した（補足図S11; 補足表S14）。 v3.0では72個（3×24）の予想ゲノムブロックのうち71個を検出し、そのほとんどが文献10で報告されたものと一致する配置になっていた（図4；補足表S15）。 v3.0では、新たにA01染色体およびA05染色体上にF (LF) およびF (MF1) という断片化したゲノムブロックが確認されたが、ref.10では観察されなかった。 また、A07番染色体上のCブロック（MF2）とA08番染色体上のBブロック（MF1）の2つの非常に小さなゲノムブロックは、v3.0では検出できませんでした（ref.10)。 また、v3.0のA08染色体上の隣接する3つのゲノムブロック（S（MF2）、T（MF2）、B（MF1））は、S/T/Bの順になっていますが、これらは文献10ではT/B/Sとして配置されていました。.
10本の染色体にあるゲノムブロックはサブゲノムLF（赤）、MF1（緑）、MF2（青）に割り当てられています。 1つのブロックの2つ以上のセグメントに小文字（a, bなど）でラベルをつけた。 B. rapaゲノムのセントロメアを黒い楕円で、古セントロメアを灰色の三角形で示した。 下向きの矢印は、1本のACK染色体に由来する他のブロックに対して反転しているGBに隣接している
また、v3.0のゲノムブロックの向きを文献10と比較しました。 A02染色体のゲノムブロックW（MF2）、A07染色体のゲノムブロックG（LF）、E（LF）は、1本のACK染色体から生じた他のブロックと比較して、逆向きであることがわかりました。 しかし、A09染色体上のゲノムブロックP（LF）とv3.0のVの3ブロックの向きは順方向であったが、ref.10ではこれらは逆方向であった。 これらの結果は、v3.0とv1.5のそれぞれの遺伝子地図でさらに支持された。
v3.0の全染色体のセントロメアの位置を正確に決定した。 セントロメアサテライトリピートCentBr, CRB, TR238, PCRBr21,22,23などの既知のセントロメアリピート配列をスクリーニングすることにより、v3.0では21の古染色体領域すべてのシグナルを確認したが、参考文献10では三つの古染色体領域を検出しなかった（図4、補足表S16）。 古染色体解析の結果、現存する10個のB. rapaのセントロメアは、すべて21個の古染色体から受け継がれたものであることがわかった。 A01、A03、A04、A05、A06、A07、A10染色体のセントロメアは、ref.10で報告されたように、対応するセントロメアを挟むゲノムブロックが同じであった（図4）。 しかし、A02染色体のセントロメアはゲノムブロックP（MF2）とV（MF1）の間に、A09染色体のセントロメアはゲノムブロックP（LF）とB（LF）の間に位置しており、これらは文献10では古染色体とみなされていた（図4）。 また、A08染色体のセントロメアは、ゲノムブロックC（MF1）とT（MF2）の間ではなく、ゲノムブロックT（MF2）とB（MF1）の間に位置していることがわかった（ref.10）。 さらに、v3.0ではセントロメア領域内に1188個の遺伝子が検出されたが、文献10では740個しか検出されなかった（補足表S17）。
v3.0でのセントロメアに関するアセンブリを評価するために、セントロメア領域の配列特性を分析した。 その結果、染色体の他の部分に比べて、セントロメア領域ではTEやセントロメア特異的な繰り返しが著しく多くマッピングされ、v3.0で注釈されたセントロメア領域では遺伝子密度や組み換え率が著しく低くなっていた（図5）。 また、ref.10で報告されたものと比較して、v3.0ではセントロメア特異的な繰り返しが多く検出された（補足表S17）。
すべてのデータはヒートマップで表現されています。 赤色は低い値、青色は高い値を示す。 a B. rapa genome v3.0の10本の染色体。 b v3.0の10本の染色体におけるTE密度（500kbスライドウィンドウ、100kbステップ） c v3.0の10本の染色体に沿ったセントロメア特異的リピートの分布（2Mbスライドウィンドウ、1Mbステップ）。 d v3.0の10本の染色体の遺伝子密度（2Mbスライドウィンドウ、1Mbステップ） e v3.0の10本の染色体に沿ったマーカー間の平均局所組み換え率（5Mbスライドウィンドウ、1Mbステップ）