NP|油菜泛基因组及生态型分化

1. 简介

油菜(Brassica napus)是世界上重要的油料作物，起源于地中海地区，大约是在7500年前由白菜(Brassica rapa)和甘蓝(Brassica oleracea)自然杂交加倍形成的异源四倍体。与白菜以及甘蓝的大量杂交，极大地扩宽了油菜的遗传多样性。为了适应自然季节变化比如春化需求、耐冬性以及光响应开花，油菜形成了不同的生态型：冬性油菜(WORs)、半冬性油菜(SWORs)以及春性油菜(SORs)。目前油菜已经公布了4个版本的参考基因组(冬油菜：Darmor-bzh,Tapidor以及半冬油菜：ZS11，NY7)。但是由于测序方法或者测序深度的限制，已公布的基因组的准确度以及完整度都不足以支撑结构变异(SV)的研究，而结构变异往往决定了遗传多样性以及重要的农艺性状。泛基因组(pan-genome)指的是一个物种所有基因的集合，包括核心基因集(core genes)以及非必需基因集(dispensable genes)。多个高质量的参考基因组有利于油菜基因组结构研究以及探究油菜不同生态型分化的遗传基础。

2. 结果

2.1 组装了八个高质量的油菜基因组

本文挑选了4个半冬性油菜(ZS11、Gangan、Zheyou7以及Shengli)、2个冬油菜(Tapidor、Quinta)以及2个春油菜(Westar、人工合成油菜No2127)共八个在花期、株型以及种子存在显著差异的品种。

采用PacBio SMRT测序以及Hi-C，作者将960.8M(95.3%)ZS11基因组挂载到19条染色体上，同时作者利用了BAC、RNA-seq、二代测序回帖以及遗传图谱验证了ZS11组装的准确度以及完整度。其它七个基因组则是根据与ZS11的共线性挂载到19条染色体上，为了验证这种挂载是可行的，作者分别对No2127进行了Hi-C以及BioNano光学图谱验证，对Westar进行了BioNano验证，发现这种基于共线性的挂载准确度是挺高的，至少这两种挂载序列一致性挺高的。当然这种挂载肯定对后续的基因组结构变异研究存在很大的影响。

作者在八个基因组中鉴定到了56.8~58.2%的重复序列，高于2014年Science发表的Darmor-bzh(34.8%)，说明此次组装的完整度是远远高于Darmor-bzh.

2.2 油菜与其他芸薹属的系统发生分析

作者利用1235个单拷贝直系同源基因，以拟南芥为外群进行系统发生分析，发现人工合成的油菜No2127的两个亚基因组明显更靠近油菜的两个祖先种基因组：白菜与甘蓝基因组。同义突变替换率(Ks)分析发现，油菜大约是在10000年前与白菜甘蓝分化开来，而白菜与甘蓝大约是在300万年前分化的。芸薹属作物大约在11个百万年前进行了一次全基因组复制事件。为了研究油菜A亚基因组的起源，作者利用210份油菜、组装的八个基因组以及199份白菜的SNP构建了进化树，发现所有的A基因组都起源于芜菁甘蓝(turnip)，PCA进一步验证了该结论。

a, Phylogenetic relationship of nine B. napus genomes and their diploid progenitors, B. rapa and B. oleracea. The phylogenetic tree is constructed on the basis of 1,235 conserved genes. The values on the branch are the substitutions between species and the nearest ancestor. WGT, whole genome triplication. b, A neighbour-joining tree of 210 B. napus accessions, eight assembled accessions and 199 B. rapa accessions. Each assembled accession was represented by a pentagram (left to right: Westar, Quinta, Tapidor, Shengli, Zheyou, Gangan, ZS11 and No2127). The layer rings indicate the group name of each clade. c, PCA plot of B. napus (n = 210) and B. rapa (n = 199) accessions. d, PCA plot of B. napus (n = 210) and B. oleracea (n = 119) accessions.

2.3 八个基因组之间的大片段结构变异

其它七个基因组与ZS11进行比对，76-84%的ZS11基因组与其它七个基因组之间存在一对一的同线性，同时共鉴定了1.87-3.93✖️10⁶高质量的SNP以及0.98-1.48✖️10⁶的InDels。

结构变异方面，鉴定到7.5-15.6Mb的倒置，其中有40个大片段(>50kb)的倒置；39.7-49.1Mb的易位，大部分是染色体间的易位，而且大部分易位是由同源交换(HE)产生的。研究发现共有16720-34158个区域，77.2-149.6Mb的序列(PAVs)是ZS11缺失的，这些区域含有2619-4810个基因，GO富集分析显示这些基因主要富集在防卫应答、信号转导以及逆境应答等通路上。

2.4 油菜泛基因组以及基因索引

作者将相对于其它七个基因组ZS11缺少的部分序列添加到ZS11基因组里，构建了一个~1.8Gb以及152185个基因的泛参考基因组。下一步 作者将九个基因组(加上Darmor-bzh)所有基因分成了四类：

• 泛基因集：九个基因组所有基因
• 核心基因集：至少在7个基因组中存在拷贝的基因
• 非必需基因集：在2-6个基因组中存在拷贝的基因
• 特异基因集：只存在一个基因组中的基因

对这四类基因集进行富集分析发现，核心基因集主要富集在一些最基本的功能比如生物合成调节、主要代谢产物的合成与运输等通路，而且86%的核心基因集基因都有已知的功能结构域；特异基因集主要富集在应激反应、逆境应答以及蛋白磷酸化等通路。针对九个基因组的基因，构建基因索引，文章还特异构建了一个 数据库网站用来浏览查阅八个基因组的所有信息。

2.5 SNP-based GWAS与PAV-based GWAS：角果长，粒重以及花期

已经有很多研究表明相比SNP，SV对性状影响更大，基于SNP的GWAS无法检测到基因组缺失对性状的影响。前人研究已经表明PAVs会影响油菜一些十分重要的性状，为了研究这种影响，作者基于3971412个SNP以及27216个PAV对角果长(SL)、粒重(SW)以及花期进行了GWAS。

基于SNP的GWAS鉴定到了75个SNP(角果长)以及38个SNP(粒重)。虽然这些SNP落在已报道的QTL区域，但是没有关联到靶基因BnaA9.CYP78A9(影响角果长以及粒重的基因)。令人振奋的是基于PAV的GWAS鉴定到了 BnaA9.CYP78A9 上游区域存在一个3.9kb的CACTA类转座子插入，而且群体中含有此转座子插入的品种表现出更长的角果以及更重的粒重，前人研究已经表明该插入的转座子作为一个增强子促进 BnaA9.CYP78A9 在角果的优先表达，进而产生更长的角果和更大的种子。在八个基因组中，Gangan、ZS11含有该转座子插入，表现出更长的角果以及更大的种子。

a, Manhattan plots of SNP-GWAS and PAV-GWAS for silique length. b, A 3.6-kb CACTA-like insertion as lead PAV of BnaA09.CYP78A9 promoter region. c, The silique length in lines with different CYP78A9 alleles. For a and b, the GWAS (-lmmm 1: Wald test) was performed with 3,971,412 SNPs or 27,216 PAVs in the BN-NAM population containing 2,141 RILs. d, Thousand-seed weight in lines with different CYP78A9 alleles. For c and d, P values were determined using two-tailed Student’s t-tests. The middle bars represent the median while the bottom and top of each box represent the 25th and 75th percentiles, respectively. The whiskers extend to 1.5 times the interquartile range. Alt, alternative; Ref, reference. e,f, Phenotype data of silique length in eight B. napus accessions. g,h, Phenotype data of seed weight in eight assembled B. napus accessions. For e and g, experiments were repeated five times with similar results. For f and h, data are mean ± s.d. of eight and five biological replicates, respectively.

花期是一个非常复杂的数量性状，为了研究花期，作者将油菜群体种植在八个环境中(2个春环境，6个冬环境)并统计花期。基于SNP的GWAS鉴定到63个SNP(春环境)以及79个SNP(冬环境)，大部分SNP是相同的，通过这些SNP，大部分调控开花的基因比如3个FLOWERING LOCUS C(FLC),2个FLOWERING LOCUS T (FT), 2个 VERNALIZATION INSENSITIVE 3 (VIN3) 以及2个CRYPTOCHROME 2 (CRY2) 都被鉴定到。基于PAV的GWAS直接关联到两个 FLC 基因: BnaA02.FLC、BnaA10.FLC ，BnaA02.FLC 的第六个外显子存在hAT(824bp)插入，群体中含有824bp的hAT的品种更早开花(春冬环境皆是)，BnaA10.FLC启动子区域存在一个4421bp的hAT插入，群体中含有该转座子插入的品种更晚开花。我们实验室曾经鉴定到了这个基因，我们当初对该基因启动子进行PCR的时候发现有些品种PCR不出来，当时我们就知道该基因上游存在结构变异，而且不是小片段的结构变异，此次PAV-GWAS验证了我们的想法。

a,b, Manhattan plots for flowering time analysed by SNP-GWAS in winter and spring environments, respectively. The gray dashed lines indicate the significance threshold. The BLUP values of the days from sowing to flowering (DTF) in the winter and spring environments were used to represent the flowering time for SNP-GWAS. The triangles and arrows denote the main candidate genes surrounding the strong peaks. c,d, Manhattan plots for flowering time analysed by PAV-GWAS in winter and spring environments, respectively. The BLUP values of DTFs in winter and spring environments were used to represent the flowering time for PAV-GWAS. The gray dashed lines indicate the significance threshold. e–g, Local Manhattan plots, gene positions and LD heatmaps show the regions surrounding the strong peaks of the candidate genes (BnaA02.FLC, BnaA10.FLC and BnaC02.FLC) identified by SNP-GWAS. h, An 824-bp hAT insertion in the last exon of BnaA02.FLC was identified as the lead PAV by PAV-GWAS. For a–h, the GWAS (-lmmm 1: Wald test) was performed with 3,971,412 SNPs or 27,216 PAVs in the BN-NAM population containing 2,141 RILs. i,j, Flowering time of lines with different BnaA02.FLC alleles in spring (i) and winter (j), respectively. P values were determined using two-tailed Student’s t-tests. The middle bars represent the median, while the bottom and top of each box represent the 25th and 75th percentiles, respectively. The whiskers extend to 1.5 times the interquartile range.

2.6 FLC基因在不同生态型油菜分化分化中的作用

FLC被报道是调控推迟开花的一个关键转录因子，它会抑制FT, SUPPRESSOR OF OVEREXPRESSION OF CO1 (SOC1)以及FD的表达。油菜含有多拷贝FLC基因以及一些跟开花相关的同源基因。作者在八个基因组以及白菜甘蓝中比较了BnaA02.FLC, BnaA10.FLC and BnaC02.FLC三个基因，发现BnaA10.FLC的编码区以及启动子区域存在四个转座子，其中一个5565-bp LINE转座子只存在于两个春油菜(No2127、Westar)的第一个外显子中；一个621-bp的MITE转座子只存在于两个冬油菜(Quinta、Tapidor)的启动子区域，这个转座子曾被报道过是专属于冬油菜的；一个4421-bp的hAT转座子只存在于四个半冬油菜(ZS11,Zheyou7,Gangan,Shengli)的启动子区域，另外还有一个1656-bp的LTR转座子存在于两个半冬油菜(Zheyou7,Shengli)的启动子区域。 进一步在210份(141份有春冬半冬分型)油菜(7X)验证这些转座子，作者选择了5个显著的序列区域进行验证，发现在这210份油菜中，34份含有MITE插入，60份含有hAT插入，129份含有LINE插入，尤其是所有的冬油菜都含有MITE插入，85%(22/26)的春油菜含有LINE插入，81%(80/99)的半冬油菜含有hAT插入。这些结果暗示着BnaA10.FLC基因中的转座子插入与否跟油菜不同生态型分化是存在十分显著的关联。在此基础上作者设计了转座子特异性引物进行PCR扩增，在BnaA10.FLC上下游5kb区域基于SNP以及转座子进行分型分析发现，转座子的插入与否与油菜生态型分化的关联是远远大于SNP的作用，因此BnaA10.FLC区域的转座子插入可以作为一种油菜生态型分型的依据之一。作者进一步在199份白菜以及192份油用卷心菜中检测BnaA10.FLC基因中的转座子插入情况，显示295份(75.4%)含有4421-bp的hAT插入以及25份(6.4%)含有MITE插入，只有18.2%的品种不含有上述任意一种转座子插入，说明这两种转座子的插入时间是早于油菜的形成时间的，春油菜特异性的LINE转座子并未在这些白菜中检测到，暗示着LINE的插入可能是在油菜的驯化选择期间发生的。 BnaA02.FLC基因第六个外显子中824-bp的hAT插入只在两个春材料中检测到。BnaC02.FLC基因在ZS11、Quinta以及Westar中存在，但是在Shengli、Zheyou7以及Tapidor中却是由于同源交换HE，该基因被BnaA02.FLC替代，至于在Gangan以及No2127中BnaC02.FLC基因完全被删除，未曾检测到。因此，Tapidor、Zheyou7以及Shengli油菜含有两个BnaA02.FLC拷贝，ZS11, Westar以及Quinta含有一个BnaA02.FLC拷贝和一个BnaC02.FLC拷贝，Gangan中只有一个BnaA02.FLC拷贝，No2127没有BnaC02.FLC拷贝。 油菜晚花以及春化反应跟BnaFLC的表达量息息相关，为了研究BnaA10.FLC,BnaA02.FLC以及BnaC02.FLC基因在不同生态型中的表达情况，从四叶期开始，每间隔一个月取一次样，在春化前以及春化后研究FLC以及FT的表达情况，发现春化前BnaA02.FLC与BnaA10.FLC在不同生态型间的表达差异十分显著，这也与由于PAV导致两基因在不同生态型间的拷贝数不同相关；春油菜(Westar，No2127)中的BnaA02.FLC基因在任何时期都不表达，BnaC02.FLC基因只在Quinta, ZS11以及Westar中表达，有研究发现BnaA02.FLC比BnaC02.FLC具有更强的抑制开花作用，虽然Tapidor与Quinta都含有BnaA10.FLC，但是Tapidor有两个BnaA02.FLC拷贝，Quinta只有一个BnaA02.FLC拷贝以及一个BnaC02.FLC拷贝，所以Tapidor与Quinta的花期存在差异。这三个FLC基因的表达量在春化之后都显著降低，同时FT基因((BnaA02.FT, BnaA07.FT以及BnaC06.FT)的表达量显著上升。FLC与FT基因的表达此消彼长，不同生态型的FLC与FT基因的拷贝数不同，表达量不同，进而导致花期差异。

a, Insertions of four transposable elements around BnaA10.FLC in different ecotypes. b, Genotyping BnaA10.FLC in 141 B. napus accessions. The left were ecotypes of B. napus accessions. The middle is the read coverage of resequencing data in 15 representative sites, with Tapidor A10: 22,661,433–22,661,437; Westar A10: 23,731,730–23,731,734 and ZS11 A10: 23,942,298–23,942,302. The right is the PCR results statistics of three insertions. c, The haplotypes of six SNPs and three TEs around BnaA10.FLC in 141 B. napus accessions. d,e, SVs in BnaA02.FLC (d) and BnaC02.FLC (e) in eight accessions. f, The expression levels of BnaA02.FLC, BnaA10.FLC and BnaC02.FLC in plants before and after vernalization (T0 and T4) based on the number of fragments per kilobase of the exon model per million mapped reads (FPKM). ***P < 0.001 (two-tailed Student’s t-test). Error bars indicate the mean ± s.d. (n = 2). ns, not spring; w, winter; s, spring; sw, semi-winter. g, The relationship between the accumulated days with low temperature, the cumulative expression levels of FLCs, FTs and flowering time in the 2018–2019 growing season in Wuhan. Three expressed FT genes were considered (average FPKM ≥ 1). Accumulated low-temperature curves indicated that the end of vernalization was in T2–T3 for SORs and T3–T4 for WORs. h, The cumulative expression levels of three FLC genes and the flowering time characterization of eight assembled B. napus accessions. Stacked histogram showed FLCs expression in T0–T3. These plants were transplanted from the field to the pot at 106 d after sowing. The standard deviation and average of flowering time were counted from 14–21 lines.

3. 主要结论

• 文章组装了八个高质量基因组，并在此基础上鉴定了大量的结构变异(SV)以及PAV，进而构建了油菜泛基因组。

• 基于PAVs的GWAS鉴定到了以前未曾鉴定到的转座子插入引起的性状差异，不同生态型之间由于PAV导致不同性状相关基因的差异表达

• 最后作者着重研究了FLC基因中的三种类型的SV，发现BnaA10.FLC, BnaA02.FLC以及BnaC02.FLC都参与了春化作用，但是BnaA10.FLC起着决定性作用。由于春油菜中BnaA10.FLC第一个外显子的LINE插入，导致春油菜中的该基因功能缺失，春油菜需要更弱或者不需要春化就能开花；冬油菜中BnaA10.FLC基因启动子区域的MITE插入上调了该基因的表达，所以冬油菜需要更强的春化作用才能开花；半冬性油菜需要的春化强度介于春油菜、冬油菜之间，这主要是由于BnaA10.FLC基因启动子区域的hAT插入导致的。

这篇文章公布的八个参考基因组质量都挺高，我个人认为有点瑕疵的是其它几个基因组是根据与ZS11的共线性挂载到染色体的，这肯定会影响到后面结构变异的分析，所以我们可以看到后面文章着重研究转座子对基因表达的影响，但是开花部分转座子的研究十分精彩，很大程度上补缺了我们实验室上篇文章的研究空白。