田志濤 ，趙永國 ，Lenka Havlickova，He Zhesi，Andrea L Harper，Ian Bancroft，鄒錫玲，張學(xué)昆，陸光遠(yuǎn)

（1中國農(nóng)業(yè)科學(xué)院油料作物研究所/農(nóng)業(yè)部油料作物生物學(xué)重點(diǎn)實(shí)驗(yàn)室，中國武漢430062；2CNAP，University of York，York YO10 5DD，UK）

0 引言

【研究意義】甘藍(lán)型油菜是世界上僅次于大豆的第二大油料作物，同時也是中國最主要的冬季油料作物，其產(chǎn)油量占國產(chǎn)油料作物產(chǎn)油量的45%以上[1-2]。油菜籽榨油后剩余的餅粕含有豐富的蛋白質(zhì)，是發(fā)展畜禽漁業(yè)的良好蛋白飼料源。然而，傳統(tǒng)油菜的種子中含有大量的硫苷，影響動物的適口性，并不適合于直接作飼料[3]。硫苷雖然本身無毒，但在芥子酶的作用下可降解為具有毒性的異硫氰酸等物質(zhì)，影響油菜餅粕作為飼料的質(zhì)量[4]。因此，降低種子硫苷含量成為世界各國油菜品質(zhì)育種的主要目標(biāo)。通過多年的努力，在20世紀(jì)90年代，油菜主要生產(chǎn)國基本普及了雙低（低芥酸、低硫苷）品種，種子中的硫苷含量被控制在30 μmol·g-1以下，產(chǎn)生了顯著的社會經(jīng)濟(jì)效益。然而，長期推廣應(yīng)用雙低油菜也出現(xiàn)了一些嚴(yán)重的問題，與傳統(tǒng)的雙高油菜相比，雙低油菜的抗病性、抗蟲性和抗逆性普遍下降，容易發(fā)生災(zāi)害，造成產(chǎn)量損失[5-6]。究其原因，是由于油菜種子中積累的硫苷主要是在葉片中合成，種子中的硫苷含量與葉片、莖桿等組織的硫苷含量高度相關(guān)，而雙低油菜選育過程并沒有打破這種相關(guān)性，在降低種子硫苷含量的同時，葉片等組織的硫苷含量也同步降低，從而導(dǎo)致油菜整體的抗病性抗逆性下降[7]。為此，有學(xué)者提出了新的品質(zhì)育種策略，即在育種過程中僅僅降低種子中的硫苷含量而保持其他組織部位較高的硫苷含量，以便保留其優(yōu)良的抗病性和抗逆性。因此，開展油菜硫苷代謝及轉(zhuǎn)運(yùn)相關(guān)遺傳因子的研究對于新時期的油菜品質(zhì)改良具有重要意義?！厩叭搜芯窟M(jìn)展】擬南芥中關(guān)于硫苷合成通路上的基因已基本明確，主要涉及到氨基酸鏈的延伸、硫苷核心結(jié)構(gòu)的生成以及側(cè)鏈的修飾[8]。硫苷在植物中的積累是動態(tài)的，其生成部位與積累部位并不相同，通常認(rèn)為，硫苷主要是在葉片中合成，后期通過莖稈轉(zhuǎn)運(yùn)到種子中。在角果成熟期，特殊的轉(zhuǎn)運(yùn)蛋白促使硫苷從角果皮轉(zhuǎn)運(yùn)到胚胎中[9]。有研究表明，擬南芥gtr1和gtr2雙突變體的種子硫苷含量極低（常規(guī)方法無法檢出），而葉片的硫苷含量則上升了10倍。這兩個基因都編碼細(xì)胞膜上的硫苷轉(zhuǎn)運(yùn)蛋白，負(fù)責(zé)將硫苷從質(zhì)外體向韌皮部轉(zhuǎn)運(yùn)[10]。在甘藍(lán)型油菜和白菜中，目前已通過傳統(tǒng)QTL分析和全基因組關(guān)聯(lián)分析發(fā)現(xiàn)了一些硫苷含量相關(guān)的QTL區(qū)間和候選基因[11-15]。然而，這些位點(diǎn)都只是通過對一個時期（成熟種子）的硫苷含量進(jìn)行分析得到，因而不能全面、準(zhǔn)確了解種子硫苷積累的遺傳機(jī)理[16-18]?！颈狙芯壳腥朦c(diǎn)】全基因組關(guān)聯(lián)分析（GWAS）是剖析復(fù)雜性狀遺傳機(jī)制的有效方法，其主要原理是充分利用歷史重組事件形成的連鎖不平衡性來推斷控制群體表型變異的遺傳位點(diǎn)[19]。因此，關(guān)聯(lián)分析的突出優(yōu)點(diǎn)是具有更加精確的定位效果和能檢測到更多的有效變異位點(diǎn)[20-21]。隨著二代測序技術(shù)的迅猛發(fā)展和測序成本的降低，基于轉(zhuǎn)錄組學(xué)的關(guān)聯(lián)分析方法應(yīng)運(yùn)而生，并開始在油菜等多倍體作物中得到應(yīng)用[22]。轉(zhuǎn)錄組關(guān)聯(lián)分析除了使用SNP標(biāo)記外，更重要的是還能使用基因表達(dá)量標(biāo)記（gene expression marker，GEM），從而檢測到更多的關(guān)聯(lián)區(qū)段，具有明顯的優(yōu)越性[23]。因此，利用該技術(shù)對動態(tài)變化中的硫苷進(jìn)行分析，有望發(fā)現(xiàn)新的硫苷代謝及轉(zhuǎn)運(yùn)相關(guān)基因。【擬解決的關(guān)鍵問題】本研究采用 113份甘藍(lán)型油菜品種構(gòu)建關(guān)聯(lián)群體，通過轉(zhuǎn)錄組測序數(shù)據(jù)開發(fā)出SNP標(biāo)記和GEM標(biāo)記，對油菜成熟種子的硫苷含量以及授粉后25 d的角果皮和種子中的硫苷含量進(jìn)行轉(zhuǎn)錄組關(guān)聯(lián)分析，以篩選與硫苷代謝及轉(zhuǎn)運(yùn)相關(guān)的候選基因，來揭示油菜硫苷變異的遺傳基礎(chǔ)，以及用于分子輔助育種來改良油菜品質(zhì)。

1 材料與方法

1.1 材料和樣品采集

本試驗(yàn)采用的甘藍(lán)型油菜材料通過項目合作從英國引進(jìn)，主要是歐洲不同時期選育的冬性品種，其中59%的材料屬于低硫苷品種（＜30 μmol·g-1）（表 1）。所有材料均經(jīng)過DH培養(yǎng)純化，一致性好，表型穩(wěn)定，在國內(nèi)已經(jīng)繁殖5代，生育期基本一致，能夠適應(yīng)長江流域氣候。

試驗(yàn)材料于2015年9月至2016年5月在中國農(nóng)科院油料所武昌基地種植，每份材料種植2行區(qū)（行長2.2 m，行距0.3 m，株距0.2 m），設(shè)置3次重復(fù)（均取樣），隨機(jī)區(qū)組排列。進(jìn)入盛花期，在每個小區(qū)選取3株長勢健壯的單株，在主花序上用棉線標(biāo)記當(dāng)天開放的花朵，然后套袋自交，謝花后及時去除袋子，讓其繼續(xù)正常生長發(fā)育。分別在授粉后的15 d（15 DAP）和25 d（25 DAP）采集鮮嫩角果，小區(qū)內(nèi)單株間混樣，立即投入液氮速凍，真空冷凍干燥后將角果皮和幼嫩種子手工剝離，4℃密封保存?zhèn)溆?。成熟時，在相同的單株上收獲種子。

1.2 硫苷檢測方法

對每次重復(fù)的混合樣品進(jìn)行硫苷檢測，最后以3次重復(fù)的平均值作為最終硫苷含量。依照[24-25]的方法提取、測定硫苷。各稱取0.1 g自然干燥的成熟油菜種子、低溫冷凍干燥的幼嫩角果皮以及幼嫩種子。將樣品加入到含有300 μL 70%甲醇溶液以及一個直徑4 mm鋼珠的2 mL離心管中。用高速研磨機(jī)將樣品研磨充分，再加入 700 μL的 70%甲醇溶液以及100 μL 濃度為 5 μmol·L-1的丙烯基硫苷標(biāo)準(zhǔn)液。超聲波提取10 min后，5 000×g離心取上清。用1 mL 70%甲醇再提取一次。將100 μL混合后的提取液加到含甲酸咪唑處理過的葡聚糖凝膠的 96孔水系過濾板中，分別用300 μL濃度為1 mol·L-1的乙酸鈉溶液沖洗兩次，加入硫酸酯酶，35℃酶解過夜。以超純水洗滌過濾板，收集濾液加載到高效液相色譜儀（安捷倫1200，美國）中，色譜條件參照國際硫苷檢測標(biāo)準(zhǔn)[25]。

1.3 轉(zhuǎn)錄組測序及SNP和GEM的檢測

RNA提取、測序、標(biāo)記開發(fā)試驗(yàn)在英國約克大學(xué)進(jìn)行。用E.Z.N.A植物RNA試劑盒（Omega Bio-tek,400 Pinnacle Way, Ste 450, Norcross, GA 30071）提取三葉期油菜幼苗的葉片RNA，具體提取方法參照試劑盒提供的說明書進(jìn)行。用焦炭酸二乙酯處理過的純水溶解提取的RNA樣品。取1 μL RNA樣品在安捷倫RNA 6000 NanoLabChip上檢測RNA的濃度和質(zhì)量。

符合質(zhì)量要求的RNA樣品在Illumine HiSeq2500平臺測序，后續(xù)的數(shù)據(jù)分析和標(biāo)記開發(fā)按照筆者課題組前期建立的方法進(jìn)行[23,26-27]。將每一個樣品的RNA-seq數(shù)據(jù)錨定到最近建立的油菜A和C泛轉(zhuǎn)錄組參考序列（pan-transcriptome resources for the Brassica A and C genomes）上[28]，通過序列比對和meta分析鑒定開發(fā)出SNP標(biāo)記。SNP標(biāo)記的質(zhì)量控制標(biāo)準(zhǔn)是：去除測序深度＜10，堿基測序質(zhì)量＜Q20，數(shù)據(jù)缺失＞0.25，等位基因數(shù)＞3的SNP位點(diǎn)。

表1 本試驗(yàn)用到的油菜材料Table 1 Brassica napus varieties used in this study

續(xù)表1 Continued table 1

續(xù)表1 Continued table 1

GEM標(biāo)記的開發(fā)方法：對油菜泛轉(zhuǎn)錄組參考序列中的116 098個CDS（coding DNA sequence）進(jìn)行數(shù)量化和均一化，計算各個CDS在關(guān)聯(lián)群體不同品種中的相對表達(dá)量（用RPKM值表示），共獲得116 098個GEM標(biāo)記。其中，53 889個標(biāo)記表達(dá)顯著（平均RPKM值＞0.4）。

1.4 SNP關(guān)聯(lián)分析

首先在PSIKO軟件包中利用所有的SNP標(biāo)記計算獲得Q矩陣[29]，然后將Q矩陣、SNP數(shù)據(jù)和表型數(shù)據(jù)導(dǎo)入到TASSEL 3.0中[30]。其次，在TASSEL中去掉最小等位頻率小于0.01的SNP數(shù)據(jù)，計算能估計品種間親緣關(guān)系的K矩陣[27]。最后，將所有的數(shù)據(jù)導(dǎo)入到一個混合線性模型（MLM）中，計算單個SNP標(biāo)記的顯著性P值和效應(yīng)值，再用R包畫出manhattan散點(diǎn)圖[22]。

1.5 GEM關(guān)聯(lián)分析

首先去除表達(dá)量（RPKM）小于0.4的GEM標(biāo)記，然后再用 GAPIT R包中的線性回歸方法計算基因表達(dá)量和GS含量的關(guān)系。對于每一個GEM標(biāo)記（即unigene），以RPKM值為自變量，GS為因變量進(jìn)行回歸分析，輸出每一個unigene的R2和顯著性P值。以每個標(biāo)記的P值取對數(shù)后的負(fù)值為縱坐標(biāo)，以標(biāo)記在染色體上的物理位置為橫坐標(biāo)，繪制曼哈頓散點(diǎn)圖，方法同上。

關(guān)聯(lián)分析顯著性采用Bonferroni閾值為標(biāo)準(zhǔn)，計算公式為：P=0.05/標(biāo)記數(shù)量。同時還采用 5%錯誤發(fā)現(xiàn)率（FDR）為標(biāo)準(zhǔn)，其含義為在最顯著關(guān)聯(lián)的標(biāo)記中有5%被認(rèn)為是假陽性。

2 結(jié)果

2.1 不同組織硫苷含量的動態(tài)變化

在關(guān)聯(lián)分析群體中，采集15 DAP、25 DAP的幼嫩角果以及自然成熟的種子，并將角果分解為角果皮和幼嫩種子兩部分后分別測定其硫苷含量，結(jié)果見表2。授粉后15 d，種子處于發(fā)育早期，干物質(zhì)較少，種子硫苷含量平均值僅為7.58 μmol·g-1，但變幅較大（變異系數(shù)為50%），從最低的1.69 μmol·g-1到最高的20.45 μmol·g-1。同樣，同一時期中角果皮的硫苷含量較低，平均值僅有4.8 μmol·g-1，但變異系數(shù)更大，達(dá)到95%。由此可見，在授粉15 d內(nèi)硫苷并沒有在角果皮及種子中大量積累。

授粉25 d后，種子發(fā)育處于中間階段，已經(jīng)積累了較多的干物質(zhì)。測定結(jié)果表明，種子和角果皮中的硫苷含量出現(xiàn)明顯上升。其中，種子的平均硫苷含量為 12.37 μmol·g-1，比 15 DAP 種子的硫苷含量升高63%，而且基因型之間的變化大（變異系數(shù)為148%），個別基因型（Norin）的含量高達(dá) 147.21 μmol·g-1，是15 DAP時期含量的8.3倍。角果皮的硫苷含量平均值為19.45 μmol·g-1，是同時期種子含量的1.6倍，也比15 DAP的角果皮含量高3.0倍，基因型之間的差異也很大。

種子成熟后（約授粉后50 d），硫苷含量總體上進(jìn)一步升高，平均為41.18 μmol·g-1，是25 DAP種子硫苷含量的3.3倍，是15 DAP種子硫苷含量的5.4 倍；變異范圍在 8.87—111.83 μmol·g-1。群體中的硫苷含量最高值也下降到 111.83 μmol·g-1，不及25 DAP種子的最高值。由此可見，硫苷的合成、運(yùn)輸和積累是一個十分復(fù)雜的動態(tài)變化過程。

相關(guān)性分析結(jié)果表明，25 DAP種子硫苷含量與同時期角果皮硫苷含量的相關(guān)系數(shù)為0.34（P＜0.05），25 DAP種子硫苷含量與成熟種子硫苷含量的相關(guān)系數(shù)為 0.26（P＜0.05），成熟種子硫苷含量與 25 DAP角果皮硫苷含量的相關(guān)系數(shù)為0.53（P＜0.01）（表3）。由此可見，油菜成熟種子中的硫苷含量明顯與前期角果皮的硫苷含量呈正相關(guān)關(guān)系，表明角果皮中的硫苷可能在后期轉(zhuǎn)運(yùn)到種子中，最后促使油菜種子中的硫苷含量大量增加。

表3 油菜不同時期、組織中的硫苷含量的相關(guān)性Table 3 The correlation of GS content in different tissues at different developing stages

25 DAP種子和25 DAP角果皮硫苷含量及成熟種子硫苷含量的變異系數(shù)都非常大，且呈連續(xù)性分布（圖1）。這表明這3個性狀受多基因控制，屬于數(shù)量性狀遺傳。因此，適合對其進(jìn)行關(guān)聯(lián)作圖分析。

圖1 硫苷含量的頻次分布圖Fig. 1 Frequency distribution of GS content

2.2 SNP關(guān)聯(lián)分析

去除次等位頻率（MAF）小于0.01的SNP標(biāo)記，得到 256 397個高質(zhì)量的 SNP標(biāo)記用于關(guān)聯(lián)分析。以校正后的Bonferroni值（P=0.05/256397=1.95×10-7）為顯著性標(biāo)準(zhǔn)，取負(fù)對數(shù)后的近似值為6.71。

對于25 DAP種子的硫苷含量，利用SNP標(biāo)記進(jìn)行關(guān)聯(lián)分析，共篩選到158個顯著關(guān)聯(lián)的SNP位點(diǎn)（圖2-a），這些SNP位點(diǎn)來源于105個CDS。除C4染色體外，其余染色體均有顯著關(guān)聯(lián)的SNP位點(diǎn)分布，其中分布較多的染色體是 A3、C1和 C3。然而，這些SNP位點(diǎn)并沒有形成明顯的關(guān)聯(lián)峰。

圖2 硫苷含量的SNP關(guān)聯(lián)分析Fig. 2 SNP association mapping of glucosinolate contents

對成熟種子的硫苷含量進(jìn)行關(guān)聯(lián)分析，總共檢測到167個顯著的SNP位點(diǎn)，來源于57個CDS（部分CDS檢測到2—4個位點(diǎn)）。這些SNP位點(diǎn)在A2、A9、C2、C7和C9染色體上形成5個非常明顯的關(guān)聯(lián)峰（圖2-b）。其中，A9染色體上標(biāo)記數(shù)量最多，達(dá)到58個，分布在19個CDS上；位于A2染色體上的標(biāo)記Cab021711.1:1953:A顯著性最高（-log10P=10.7），可解釋表型變異的68%（表4）。相較于25 DAP種子硫苷的SNP關(guān)聯(lián)分析，成熟種子硫苷的SNP關(guān)聯(lián)分析檢測到的關(guān)聯(lián)位點(diǎn)更加集中。

為了探討角果皮在硫苷合成和轉(zhuǎn)運(yùn)中的作用，對授粉后25天的角果皮硫苷含量也進(jìn)行SNP關(guān)聯(lián)分析（圖2-c）。結(jié)果表明，僅在A8、A9和C3染色體上各檢測到1個顯著關(guān)聯(lián)位點(diǎn)，遠(yuǎn)遠(yuǎn)少于種子的檢測結(jié)果。

表4 成熟種子硫苷含量SNP關(guān)聯(lián)分析結(jié)果Table 4 SNP association peak for glucosinolate content in mature seeds

2.3 GEM關(guān)聯(lián)分析

去除表達(dá)量（RPKM）小于0.4的 GEM 標(biāo)記,得到53 889個高質(zhì)量GEM標(biāo)記，其中，25 834個位于A基因組，28 055個位于C基因組。以硫苷含量為因變量，GEM標(biāo)記為自變量，利用混合線性模型進(jìn)行關(guān)聯(lián)分析。以校正后的 Bonferroni值（P=0.05/53889=9.28×10-7）為顯著性標(biāo)準(zhǔn)，取負(fù)對數(shù)后的近似值為6.03。另外，在控制假陽性率的情況下，利用多重假設(shè)檢驗(yàn)FDR=0.05來檢測可能由于Bonferroni值太過嚴(yán)格而遺漏的候選基因或關(guān)聯(lián)區(qū)間。

結(jié)果顯示，與25 DAP種子硫苷含量顯著關(guān)聯(lián)的GEM標(biāo)記有16個，分布于A2、A3、A8、A9、C2、C3、C5和C9染色體上，其中分布在A8和C2染色體上的標(biāo)記較多，分別為5個和3個，其余染色體僅有 1—2個標(biāo)記（圖 3-a）。若以 FDR檢驗(yàn)為標(biāo)準(zhǔn)（-log10P=4.5），則篩選到的標(biāo)記數(shù)量達(dá)到45個。

成熟種子硫苷含量共檢測到127個顯著GEM標(biāo)記，分布于除A4和A10染色體外的其他所有17條染色體上（圖3-b），其中在A2、A9、C2和C9染色體上出現(xiàn)了4個明顯的關(guān)聯(lián)峰（表5），即GEMb-A2、GEMb-A9、GEMb-C2和 GEMb-C9，其中位于GEMb-C2上的標(biāo)記數(shù)量最多（標(biāo)記Bo2g161790.1的顯著性高達(dá)-log10P=23.42，貢獻(xiàn)率為 65%）；其次是GEMb-C9，含有 26個標(biāo)記，顯著性最高為-log10P=10.15；標(biāo)記最少的峰是GEMb-A2，僅有13個，且代表性標(biāo)記Cab021665.1的顯著性僅為-log10P=12.1。

表5 成熟種子及授粉25天角果皮硫苷含量GEM關(guān)聯(lián)分析結(jié)果Table 5 Association mapping of glucosinolate contents in mature seeds and 25 DAP silique walls with GEM

授粉25 d后角果皮硫苷含量共檢測到24個顯著性標(biāo)記，分布于 A1、A2、A8、A9、C2、C4和 C8染色體上（圖3-c）。C9染色體上有一個明顯的關(guān)聯(lián)峰，但未達(dá)Bonferroni顯著性閾值，因此采用假陽性率FDR=0.05為新的閾值來篩選關(guān)聯(lián)區(qū)間，最終得到GEMc-A8、GEMc-A9、GEMc-C2和GEMc-C9四個關(guān)聯(lián)區(qū)間（表5）。

圖3 硫苷含量的GEM關(guān)聯(lián)分析Fig. 3 GEM association mapping of glucosinolate contents

2.4 不同標(biāo)記關(guān)聯(lián)分析比較

GEM關(guān)聯(lián)分析表明，25 DAP角果皮與成熟種子的硫苷含量共同檢測到 3個明顯重合的關(guān)聯(lián)區(qū)間，分別位于A9、C2和C9上（圖3和表5）。不僅在區(qū)間上存在重合，而且有22個顯著標(biāo)記在兩個性狀間也是相同的。而成熟種子與25 DAP種子的硫苷含量卻只檢測到1個共同的顯著標(biāo)記，25 DAP種子與25 DAP角果皮的硫苷則沒有檢測到共有的顯著標(biāo)記。

作為對比，在SNP關(guān)聯(lián)分析中，各個性狀之間沒有檢測到共同的標(biāo)記。但是，以顯著性SNP標(biāo)記所在的CDS為目標(biāo)，SNP標(biāo)記關(guān)聯(lián)分析共篩選出165個CDS，其中有14個同時被GEM標(biāo)記關(guān)聯(lián)分析檢測到（共篩選出的144個CDS）（表6）。

另外，本研究檢測到的主要關(guān)聯(lián)區(qū)域與前人報道結(jié)果相吻合[13-14,31-32]（表7）。比如，F(xiàn)ENG等[31]利用 DH群體定位到的油菜種子硫苷含量的 4個QTL，以及LI等[33]利用472份油菜進(jìn)行關(guān)聯(lián)分析獲得的3個關(guān)聯(lián)區(qū)段，在本研究中得到進(jìn)一步證實(shí)。以上比對結(jié)果說明，本研究所采用的技術(shù)策略是有效的，而且比前人檢測到更多顯著的關(guān)聯(lián)標(biāo)記，因而功能更加強(qiáng)大。

表6 SNP和GEM共同檢測到的CDSTable 6 Common CDS detected by SNP and GEM in association mapping

表7 不同油菜成熟種子硫苷含量定位研究的比較Table 7 The comparison of different traits on genetic analysis of GS content of mature rapeseed

2.5 候選基因的預(yù)測

以擬南芥基因功能注釋為參考，通過序列相似性比對來篩選與硫苷代謝相關(guān)的油菜同源基因。整合不同時期、不同組織器官硫苷含量檢測到的關(guān)聯(lián)區(qū)間，共得到6個關(guān)聯(lián)區(qū)間，分別位于A2、A8、A9、C2、C8和 C9染色體上。在關(guān)聯(lián)區(qū)間內(nèi)提取所有顯著的CDS，將 CDS的序列比對到擬南芥數(shù)據(jù)庫（Tair，www.arabidopsis.org），并進(jìn)行功能注釋（表 8）。在A2、A9、C2和C9染色體上共發(fā)現(xiàn)15個參與硫苷代謝的候選基因。其中A2染色體關(guān)聯(lián)區(qū)間內(nèi)所含候選基因最少，只有2個候選基因，而其他區(qū)間都含有4—5個候選基因。在4個關(guān)聯(lián)區(qū)間內(nèi)，MYB28轉(zhuǎn)錄因子蛋白被重復(fù)檢測到，而且該轉(zhuǎn)錄因子在以往關(guān)聯(lián)分析研究中也被多次檢測出與種子硫苷存在明顯的關(guān)聯(lián)性[22-23]。而其他基因則在硫元素同化、氨基酸合成、硫苷核心結(jié)構(gòu)生成以及硫苷轉(zhuǎn)運(yùn)過程中發(fā)揮重要作用。

表8 SNP和GEM關(guān)聯(lián)峰內(nèi)篩選到的候選基因Table 8 Summary of candidate genes within association peaks with SNPs and GEMs

表9 關(guān)聯(lián)峰區(qū)間外篩選到的候選基因Table 9 Summary of candidate genes out of association peak intervals

在關(guān)聯(lián)峰區(qū)間外提取兩種顯著標(biāo)記（SNP或者GEM）所在的CDS，并進(jìn)行功能注釋，篩選出與硫苷代謝相關(guān)的候選基因（表9）。以SNP為標(biāo)記共檢測到3個硫苷代謝相關(guān)候選基因。其中，BnWRKY21和BnMYB73擬南芥突變體表現(xiàn)為明顯的低硫苷水平[33]，而 BnCYP81D3則是預(yù)測為吲哚族硫苷代謝相關(guān)的基因。這3個候選基因都是在25 DAP種子的關(guān)聯(lián)分析中檢測到的，分布在A3和C1染色體上。GEM關(guān)聯(lián)分析僅在成熟種子中檢測到6個硫苷代謝相關(guān)候選基因，分布在 A3、A6、C5和 C7等多個染色體上。

3 討論

利用一個歐洲甘藍(lán)型冬油菜自然群體，首次開展了不同組織器官硫苷含量動態(tài)變化分析和轉(zhuǎn)錄組關(guān)聯(lián)分析。硫苷主要分布在成熟種子以及授粉25 d后的角果皮中，而幼嫩種子中的硫苷含量很少。相關(guān)分析表明，成熟種子中硫苷積累量與角果皮中硫苷的含量存在明顯的正相關(guān)關(guān)系。在GEM標(biāo)記的關(guān)聯(lián)分析中，25 DAP角果皮檢測到的顯著性GEM標(biāo)記大部分（92%）在成熟種子中能重復(fù)檢測到，暗示油菜中的硫苷也可能像擬南芥、蕪菁以及芥菜一樣存在由角果皮到種子的轉(zhuǎn)運(yùn)過程[10,34]。

甘藍(lán)型油菜具有A和C兩個染色體組，分別來源于白菜和甘藍(lán)。這兩套遺傳物質(zhì)非常相似，其中在結(jié)構(gòu)上存在大約15%的差異，而表達(dá)模式僅存在3%的差異[22,35]。這種染色體間高度的序列相似性嚴(yán)重影響到全基因組關(guān)聯(lián)分析的檢測效果。為此，本研究利用了一種新的關(guān)聯(lián)分析技術(shù)—轉(zhuǎn)錄組關(guān)聯(lián)分析（associative transcriptomic，AT），這種技術(shù)既可以產(chǎn)生高密度的SNP標(biāo)記，同時也能產(chǎn)生大量的GEM標(biāo)記。利用這兩種標(biāo)記來分析控制表型性狀的遺傳機(jī)制能克服全基因組關(guān)聯(lián)分析在多倍體中應(yīng)用的不足[23]。有研究表明，大量基因在幼嫩葉片中表達(dá)，葉片是營養(yǎng)生長階段硫苷合成的主要部位，在生殖生長階段再轉(zhuǎn)運(yùn)至植株的其他組織器官中[10,24,36-37]。因此，葉片上表達(dá)的基因與后期角果皮以及種子中硫苷的積累存在生物學(xué)上的聯(lián)系。在本研究中，利用油菜嫩葉提取的RNA，通過轉(zhuǎn)錄組測序技術(shù)開發(fā)了大量的SNP和 GEM標(biāo)記，在此基礎(chǔ)上對油菜兩個發(fā)育時期的種子以及角果皮硫苷含量進(jìn)行轉(zhuǎn)錄組關(guān)聯(lián)分析。結(jié)果發(fā)現(xiàn)，在成熟種子中共檢測到5個明顯的關(guān)聯(lián)峰，其位置與前人的 QTL作圖分析以及普通全基因組關(guān)聯(lián)分析結(jié)果相吻合[13-14,31-32]。這也表明，轉(zhuǎn)錄組關(guān)聯(lián)分析方法具有很好的可靠性和穩(wěn)定性。利用SNP和GEM標(biāo)記分別檢測出165和144個候選CDS，其中有14個CDS被兩種標(biāo)記共同檢測到。通過擬南芥同源基因功能注釋，發(fā)現(xiàn) GEM檢測到與硫苷直接相關(guān)候選基因的比例明顯多于SNP標(biāo)記。這進(jìn)一步說明，轉(zhuǎn)錄組數(shù)據(jù)提供額外的 GEM標(biāo)記能對全基因組關(guān)聯(lián)分析方法起到驗(yàn)證作用和補(bǔ)充效果。

基于動態(tài)表型性狀，關(guān)聯(lián)分析可以充分挖掘控制表型性狀的遺傳位點(diǎn)[16]。本研究中，在成熟種子中共掃描到5個關(guān)聯(lián)峰，在25 DAP角果皮中共掃描到4個關(guān)聯(lián)峰，而在25 DAP種子中沒有掃描到明顯的關(guān)聯(lián)峰。以上結(jié)果說明，僅僅進(jìn)行單一發(fā)育時期或單一組織器官的關(guān)聯(lián)分析可能會遺漏掉很多重要的遺傳位點(diǎn)，而這些差異的遺傳位點(diǎn)正是解釋表型動態(tài)變化的關(guān)鍵位點(diǎn)[17]。以上結(jié)果還說明，不同時期表型的變異受到不同遺傳位點(diǎn)調(diào)控[18]。利用 GEM標(biāo)記，在成熟種子與嫩角果皮中檢測到22個共同的候選基因，其中大部分候選基因與硫苷合成相關(guān)。其原因可能是，種子與角果皮中的硫苷都是在其他部位（如葉片）合成后運(yùn)輸過來的。

本研究中，成熟種子和角果皮中均檢測到3個與硫苷顯著關(guān)聯(lián)的標(biāo)記（Cab021728.1、Cab038298.3和Bo2g161590.1）。這3個標(biāo)記對應(yīng)的基因有著明顯的序列相似性，均編碼轉(zhuǎn)錄因子蛋白（Myb28），它能通過調(diào)控多個硫苷通路上的基因來影響脂肪族硫苷的合成[22-23]。在成熟種子中還檢測到多個表達(dá) MAM1蛋白的基因，其氨基酸序列與2-異丙基蘋果酸合成酶序列高度相似，在硫苷合成過程中控制脂肪族硫苷R鏈的延伸[38]。在大白菜中超表達(dá)擬南芥MAM1能導(dǎo)致3-丁烯基硫代葡萄糖苷以及 4-戊希基硫代葡萄糖苷的積累[39]。本研究還檢測到硫苷主鏈合成過程中的3個基因，CYP79F1、CYP79F2 和 UGT74B1。其中CYP79F1和CYP79F2參與脂肪族硫苷主鏈合成的第一步，這些蛋白能催化甲硫氨酸同系物轉(zhuǎn)化為相應(yīng)的醛肟[40-41]。UGT74B1則表達(dá)一種UDP依賴性葡萄糖基轉(zhuǎn)移酶。擬南芥功能互補(bǔ)試驗(yàn)表明UGT74B1控制硫苷的積累并影響植物生長素的內(nèi)穩(wěn)態(tài)[42-43]。AOP2能介導(dǎo)甲基亞砜烷基硫苷到烯烴基硫苷的轉(zhuǎn)化過程，該基因產(chǎn)物除了能參與這種轉(zhuǎn)化過程，還能促進(jìn)MYB28和MYB29這兩個控制硫苷合成的轉(zhuǎn)錄因子的表達(dá)[44-45]。雖然候選基因中還有很多并不是硫苷合成通路上的基因，但是這些基因大多屬于轉(zhuǎn)錄因子、信號響應(yīng)因子、具有氧化還原活性的酶和參與細(xì)胞過程的蛋白等。這些基因可以通過間接的方式潛在地影響硫苷的積累。例如位于A9染色體上的 APK就是通過參與硫元素的同化過程從而間接影響硫苷的合成[46]。而在A9和C9染色體的關(guān)聯(lián)區(qū)間內(nèi)找到的GTR2則表達(dá)一種硫苷轉(zhuǎn)運(yùn)蛋白，這種蛋白參與硫苷的轉(zhuǎn)運(yùn)過程，最終影響硫苷在植物種子中的積累[10,34]。這些都可能間接影響油菜硫苷在不同部位的積累。以上結(jié)果足以說明轉(zhuǎn)錄組關(guān)聯(lián)分析方法的強(qiáng)大功能。除此之外，本研究還同時檢測到多個新的候選基因，這些基因的功能還不是很明確或是還沒有證實(shí)與硫苷的積累有關(guān)，需要通過后續(xù)研究來進(jìn)一步解析其在油菜硫苷積累中的作用。

4 結(jié)論

通過對113份甘藍(lán)型油菜種子及角果皮硫苷含量的轉(zhuǎn)錄組關(guān)聯(lián)分析，獲得了328個顯著關(guān)聯(lián)的單核苷酸多態(tài)性標(biāo)記（SNP）和 144個基因表達(dá)量標(biāo)記（GEM），預(yù)測到25個與油菜硫苷代謝相關(guān)的候選基因，以及73個功能未知的新基因，這些結(jié)果可用于油菜硫苷含量的遺傳調(diào)控以及分子標(biāo)記輔助選擇研究。

[1] 張永霞, 趙鋒, 張紅玲. 中國油菜產(chǎn)業(yè)發(fā)展現(xiàn)狀問題及對策分析.世界農(nóng)業(yè), 2015, 4: 96-99, 203-204.ZHANG Y X, ZHAO F, ZHANG H L. The current situation, problem and solutions of China’s rape industry. World Agriculture, 2015, 4:96-99, 203-204. (in Chinese)

[2] 王漢中, 殷艷. 我國油料產(chǎn)業(yè)形勢分析與發(fā)展對策建議. 中國油料作物學(xué)報, 2014, 36(3): 414-421.WANG H Z, YIN Y. Analysis and strategy for oil crop industry in China. Chinese Journal of Oil Crop Sciences, 2014, 36(3): 414-421.(in Chinese)

[3] 申躍宇, 尹佩輝, 莫放. 菜籽餅(粕)的抗?fàn)I養(yǎng)因子及其對奶牛生產(chǎn)的影響. 中國草食動物, 2009, 29(2): 49-51.SHEN Y Y, YIN P H, MO F. Effect of anti-nutritional factors in rapeseed meal on dairy. China Herbivores, 2009, 29(2): 49-51. (in Chinese)

[4] TRIPATHI M K, MISHRA A S. Glucosinolates in animal nutrition: A review. Animal Feed Science and Technology, 2007, 132(1): 1-27.

[5] SCHNUG E, CEYNOWA J. Phytopathological aspects of glucosinolates in oilseed rape. Journal of Agronomy and Crop Science-Zeitschrift Fur Acker Und Pflanzenbau, 1990, 165(5): 319-328.

[6] 宋志榮, 官春云. 甘藍(lán)型油菜硫苷特性與對菌核病抗性關(guān)系. 湖南農(nóng)業(yè)大學(xué)學(xué)報(自然科學(xué)版), 2008, 34: 462-465.SONG Z R, GUAN C Y. Relationship between glucosinolate characteristics and resistance to Sclerotinia sclerotiorum in Brassica napus L. Journal of Hunan Agricultural University (Natural Science Edition), 2008, 34: 462-465. (in Chinese)

[7] LI Y C, KIDDLE G, BENNETT R, DOUGHTY K, WALLSGROVE R. Variation in the glucosinolate content of vegetative tissues of Chinese lines of Brassica napus L. Annals of Applied Biology, 1999,134(1): 131-136.

[8] HALKIER B A, DU L C. The biosynthesis of glucosinolates. Trends in Plant Science, 1997, 2(11): 425-431.

[9] DU L C, HALKIER B A. Biosynthesis of glucosinolates in the developing silique walls and seeds of Sinapis alba. Phytochemistry,1998, 48(7): 1145-1150.

[10] NOUR-ELDIN H H, ANDERSEN T G, BUROW M, MADSEN S R,JORGENSEN M E, OLSEN C E, DREYER I, HEDRICH R,GEIGER D, HALKIER B A. NRT/PTR transporters are essential for translocation of glucosinolate defence compounds to seeds. Nature,2012, 488(7412): 531-534.

[11] XU J F, LONG Y, WU J G, XU H M, ZHAO Z G, WEN J, MENG J L,SHI C H. QTL identification on two genetic systems for rapeseed glucosinolate and erucic acid contents over two seasons. Euphytica,2015, 205(3): 647-657.

[12] GAJARDO H A, WITTKOP B, SOTO-CERDA B, HIGGINS E E,PARKIN I A P, SNOWDON R J, FEDERICO M L, INIGUEZ-LUY F L. Association mapping of seed quality traits in Brassica napus L.using GWAS and candidate QTL approaches. Molecular Breeding,2015, 35(6): 143-152.

[13] LI F, CHEN B Y, XU K, WU J F, SONG W L, BANCROFT I,HARPER A L, TRICK M, LIU S Y, GAO G Z, WANG N, YAN G X,QIAO J W, LI J, LI H, XIAO X, ZHANG T Y, WU X M.Genome-wide association study dissects the genetic architecture of seed weight and seed quality in rapeseed (Brassica napus L.). DNA Research, 2014, 21(4): 355-367.

[14] LOU P, ZHAO J, HE H, HANHART C, DEL CARPIO D P,VERKERK R, CUSTERS J, KOORNNEEF M, BONNEMA G.Quantitative trait loci for glucosinolate accumulation in Brassica rapa leaves. New Phytologist, 2008, 179(4): 1017-1032.

[15] HASAN M, FRIEDT W, PONS-KUEHNEMANN J, FREITAG N M,LINK K, SNOWDON R J. Association of gene-linked SSR markers to seed glucosinolate content in oilseed rape (Brassica napus ssp. napus).Theoretical and Applied Genetics, 2008, 116(8): 1035-1049.

[16] JIANG L, LIU J, ZHU X, YE M, SUN L, LACAZE X, WU R.2HiGWAS: a unifying high-dimensional platform to infer the global genetic architecture of trait development. Briefings in Bioinformatics,2015, 16(6): 905-911.

[17] BAC-MOLENAAR J A, VREUGDENHIL D, GRANIER C,KEURENTJES J J. Genome-wide association mapping of growth dynamics detects time-specific and general quantitative trait loci.Journal of Experimental Botany, 2015, 66(18): 5567-5580.

[18] WURSCHUM T, LIU W, ALHEIT K V, TUCKER M R, GOWDA M,WEISSMANN E A, HAHN V, MAURER H P. Adult plant development in triticale (x Triticosecale Wittmack) is controlled by dynamic genetic patterns of regulation. G3 (Bethesda, Md.), 2014,4(9): 1585-1591.

[19] HALL D, TEGSTR?M C, INGVARSSON P K. Using association mapping to dissect the genetic basis of complex traits in plants.Briefings in Functional Genomics, 2010, 9(2): 157-165.

[20] YAN J, WARBURTON M, CROUCH J. Association mapping for enhancing maize (Zea mays L.) genetic improvement. Crop Science,2011, 51(2): 433-449.

[21] ZHU C, GORE M, BUCKLER E S, YU J. Status and prospects of association mapping in plants. Plant Genome, 2008, 1(1): 5-20.

[22] HARPER A L, TRICK M, HIGGINS J, FRASER F, CLISSOLD L,WELLS R, HATTORI C, WERNER P, BANCROFT I. Associative transcriptomics of traits in the polyploid crop species Brassica napus.Nature Biotechnology, 2012, 30(8): 798-802.

[23] LU G Y, HARPER A L, TRICK M, MORGAN C, FRASER F,O'NEILL C, BANCROFT I. Associative transcriptomics study dissects the genetic architecture of seed glucosinolate content in Brassica napus. DNA Research, 2014, 21(6): 613-625.

[24] KLIEBENSTEIN D J, KROYMANN J, BROWN P, FIGUTH A,PEDERSEN D, GERSHENZON J, MITCHELL-OLDS T. Genetic control of natural variation in Arabidopsis glucosinolate accumulation.Plant Physiology, 2001, 126(2): 811-825.

[25] ORIS. Determination of glucosinolates content-Part 1: Method using high performance liquid chromatography//International Standard Organization. Geneva, Switzerland, 1992: 1-9.

[26] HARPER A L, MCKINNEY L V, NIELSEN L R, HAVLICKOVA L,LI Y, TRICK M, FRASER F, WANG L H, FELLGETT A, SOLLARS E S A, JANACEK S H, DOWNIE J A, BUGGS R J A, KJAER E D,BANCROFT I. Molecular markers for tolerance of European ash(Fraxinus excelsior) to dieback disease identified using associative transcriptomics. Scientific Reports, 2016, 6: 19335.

[27] WOOD I P, PEARSON B M, GARCIA-GUTIERREZ E,HAVLICKOVA L, HE Z, HARPER A L, BANCROFT I, WALDRON K W. Carbohydrate microarrays and their use for the identification of molecular markers for plant cell wall composition. Proceedings of the National Academy of Sciences, 2017, 114(26): 6860-6865.

[28] HE Z, CHENG F, LI Y, WANG X, PARKIN I A, CHALHOUB B,LIU S, BANCROFT I. Construction of Brassica A and C genomebased ordered pan-transcriptomes for use in rapeseed genomic research. Data in Brief, 2015, 4: 357-362.

[29] POPESCU A A, HARPER A L, TRICK M, BANCROFT I, HUBER K T. A novel and fast approach for population structure inference using kernel-PCA and optimization. Genetics, 2014, 198(4): 1421-1431.

[30] BRADBURY P J, ZHANG Z, KROON D E, CASSTEVENS T M,RAMDOSS Y, BUCKLER E S. TASSEL: Software for association mapping of complex traits in diverse samples. Bioinformatics, 2007,23: 2633-2635.

[31] FENG J, LONG Y, SHI L, SHI J Q, BARKER G, MENG J L.Characterization of metabolite quantitative trait loci and metabolic networks that control glucosinolate concentration in the seeds and leaves of Brassica napus. New Phytologist, 2012, 193(1): 96-108.

[32] FU Y, LU K, QIAN L W, MEI J Q, WEI D Y, PENG X H, XU X F, LI J N, FRAUEN M, DREYER F, SNOWDON R J, QIAN W.Development of genic cleavage markers in association with seed glucosinolate content in canola. Theoretical and Applied Genetics,2015, 128(6): 1029-1037.

[33] LI B, GAUDINIER A, TANG M, TAYLOR-TEEPLES M, NHAM N T, GHAFFARI C, BENSON D S, STEINMANN M, GRAY J A,BRADY S M, KLIEBENSTEIN D J. Promoter-based integration in plant defense regulation. Plant Physiology, 2014, 166(4): 1803-1820.

[34] NOUR-ELDIN H H, MADSEN S R, ENGELEN S, JORGENSEN M E, OLSEN C E, ANDERSEN J S, SEYNNAEVE D, VERHOYE T,FULAWKA R, DENOLF P, HALKIER B A. Reduction of antinutritional glucosinolates in Brassica oilseeds by mutation of genes encoding transporters. Nature Biotechnology, 2017, 35(4):377-382.

[35] HIGGINS J, MAGUSIN A, TRICK M, FRASER F, BANCROFT I.Use of mRNA-seq to discriminate contributions to the transcriptome from the constituent genomes of the polyploid crop species Brassica napus. BMC Genomics, 2012, 13: 247-251.

[36] CHEN S X, PETERSEN B L, OLSEN C E, SCHULZ A, HALKIER B A. Long-distance phloem transport of glucosinolates in Arabidopsis.Plant Physiology, 2001, 127(1): 194-201.

[37] CHEN S, ANDREASSON E. Update on glucosinolate metabolism and transport. Plant Physiology and Biochemistry, 2001, 39(9): 743-758.

[38] KROYMANN J, TEXTOR S, TOKUHISA J G, FALK K L,BARTRAM S, GERSHENZON J, MITCHELL-OLDS T. A gene controlling variation in Arabidopsis glucosinolate composition is part of the methionine chain elongation pathway. Plant Physiology, 2001,127(3): 1077-1088.

[39] ZANG Y X, KIM J H, PARK Y D, KIM D H, HONG S B. Metabolic engineering of aliphatic glucosinolates in Chinese cabbage plants expressing Arabidopsis MAM1, CYP79F1, and CYP83A1. BMB Reports, 2008, 41(6): 472-478.

[40] MEWIS I, TOKUHISA J G, SCHULTZ J C, APPEL H M, ULRICHS C, GERSHENZON J. Gene expression and glucosinolate accumulation in Arabidopsis thaliana in response to generalist and specialist herbivores of different feeding guilds and the role of defense signaling pathways. Phytochemistry, 2006, 67(22): 2450-2462.

[41] HANSEN C H, WITTSTOCK U, OLSEN C E, HICK A J, PICKETT J A, HALKIER B A. Cytochrome P450CYP79F1 from Arabidopsis catalyzes the conversion of dihomomethionine and trihomomethionine to the corresponding aldoximes in the biosynthesis of aliphatic glucosinolates. Journal of Biological Chemistry, 2001, 276(14):11078-11085.

[42] MARROUN S, MONTAUT S, MARQUES S, LAFITE P, COADOU G, ROLLIN P, JOUSSET G, SCHULER M, TATIBOUET A,OULYADI H, DANIELLOU R. UGT74B1 from Arabidopsis thaliana as a versatile biocatalyst for the synthesis of desulfoglycosinolates.Organic & Biomolecular Chemistry, 2016, 14(26): 6252-6261.

[43] GRUBB C D, ZIPP B J, LUDWIG-MULLER J, MASUNO M N,MOLINSKI T F, ABEL S. Arabidopsis glucosyltransferase UGT74B1 functions in glucosinolate biosynthesis and auxin homeostasis. Plant Journal, 2004, 40(6): 893-908.

[44] BUROW M, ATWELL S, FRANCISCO M, KERWIN R E,HALKIER B A, KLIEBENSTEIN D J. The Glucosinolate biosynthetic gene AOP2 mediates feed-back regulation of jasmonic acid signaling in Arabidopsis. Molecular Plant, 2015, 8(8): 1201-1212.

[45] NEAL C S, FREDERICKS D P, GRIFFITHS C A, NEALE A D. The characterisation of AOP2: A gene associated with the biosynthesis of aliphatic alkenyl glucosinolates in Arabidopsis thaliana. BMC Plant Biology, 2010, 10: 170.

[46] YATUSEVICH R, MUGFORD S G, MATTHEWMAN C,GIGOLASHVILI T, FRERIGMANN H, DELANEY S, KOPRIVOVA A, FLUGGE U I, KOPRIVA S. Genes of primary sulfate assimilation are part of the glucosinolate biosynthetic network in Arabidopsis thaliana. Plant Journal, 2010, 62(1): 1-11.