甘藍(lán)型油菜種子和角果皮中硫苷含量的動(dòng)態(tài)變化及轉(zhuǎn)錄組關(guān)聯(lián)分析

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

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

0 引言

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

1 材料與方法

1.1 材料和樣品采集

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

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

1.2 硫苷檢測(cè)方法

對(duì)每次重復(fù)的混合樣品進(jìn)行硫苷檢測(cè)，最后以3次重復(fù)的平均值作為最終硫苷含量。依照[24-25]的方法提取、測(cè)定硫苷。各稱取0.1 g自然干燥的成熟油菜種子、低溫冷凍干燥的幼嫩角果皮以及幼嫩種子。將樣品加入到含有300 μL 70%甲醇溶液以及一個(gè)直徑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混合后的提取液加到含甲酸咪唑處理過(guò)的葡聚糖凝膠的 96孔水系過(guò)濾板中，分別用300 μL濃度為1 mol·L-1的乙酸鈉溶液沖洗兩次，加入硫酸酯酶，35℃酶解過(guò)夜。以超純水洗滌過(guò)濾板，收集濾液加載到高效液相色譜儀（安捷倫1200，美國(guó)）中，色譜條件參照國(guó)際硫苷檢測(cè)標(biāo)準(zhǔn)[25]。

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

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

符合質(zhì)量要求的RNA樣品在Illumine HiSeq2500平臺(tái)測(cè)序，后續(xù)的數(shù)據(jù)分析和標(biāo)記開(kāi)發(fā)按照筆者課題組前期建立的方法進(jìn)行[23,26-27]。將每一個(gè)樣品的RNA-seq數(shù)據(jù)錨定到最近建立的油菜A和C泛轉(zhuǎn)錄組參考序列（pan-transcriptome resources for the Brassica A and C genomes）上[28]，通過(guò)序列比對(duì)和meta分析鑒定開(kāi)發(fā)出SNP標(biāo)記。SNP標(biāo)記的質(zhì)量控制標(biāo)準(zhǔn)是：去除測(cè)序深度＜10，堿基測(cè)序質(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)記的開(kāi)發(fā)方法：對(duì)油菜泛轉(zhuǎn)錄組參考序列中的116 098個(gè)CDS（coding DNA sequence）進(jìn)行數(shù)量化和均一化，計(jì)算各個(gè)CDS在關(guān)聯(lián)群體不同品種中的相對(duì)表達(dá)量（用RPKM值表示），共獲得116 098個(gè)GEM標(biāo)記。其中，53 889個(gè)標(biāo)記表達(dá)顯著（平均RPKM值＞0.4）。

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

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

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

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

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

2 結(jié)果

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

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

授粉25 d后，種子發(fā)育處于中間階段，已經(jīng)積累了較多的干物質(zhì)。測(cè)定結(jié)果表明，種子和角果皮中的硫苷含量出現(xiàn)明顯上升。其中，種子的平均硫苷含量為 12.37 μmol·g-1，比 15 DAP 種子的硫苷含量升高63%，而且基因型之間的變化大（變異系數(shù)為148%），個(gè)別基因型（Norin）的含量高達(dá) 147.21 μmol·g-1，是15 DAP時(shí)期含量的8.3倍。角果皮的硫苷含量平均值為19.45 μmol·g-1，是同時(shí)期種子含量的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種子的最高值。由此可見(jiàn)，硫苷的合成、運(yùn)輸和積累是一個(gè)十分復(fù)雜的動(dòng)態(tài)變化過(guò)程。

相關(guān)性分析結(jié)果表明，25 DAP種子硫苷含量與同時(shí)期角果皮硫苷含量的相關(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）。由此可見(jiàn)，油菜成熟種子中的硫苷含量明顯與前期角果皮的硫苷含量呈正相關(guān)關(guān)系，表明角果皮中的硫苷可能在后期轉(zhuǎn)運(yùn)到種子中，最后促使油菜種子中的硫苷含量大量增加。

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

25 DAP種子和25 DAP角果皮硫苷含量及成熟種子硫苷含量的變異系數(shù)都非常大，且呈連續(xù)性分布（圖1）。這表明這3個(gè)性狀受多基因控制，屬于數(shù)量性狀遺傳。因此，適合對(duì)其進(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個(gè)高質(zhì)量的 SNP標(biāo)記用于關(guān)聯(lián)分析。以校正后的Bonferroni值（P=0.05/256397=1.95×10-7）為顯著性標(biāo)準(zhǔn)，取負(fù)對(duì)數(shù)后的近似值為6.71。

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

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

對(duì)成熟種子的硫苷含量進(jìn)行關(guān)聯(lián)分析，總共檢測(cè)到167個(gè)顯著的SNP位點(diǎn)，來(lái)源于57個(gè)CDS（部分CDS檢測(cè)到2—4個(gè)位點(diǎn)）。這些SNP位點(diǎn)在A2、A9、C2、C7和C9染色體上形成5個(gè)非常明顯的關(guān)聯(lián)峰（圖2-b）。其中，A9染色體上標(biāo)記數(shù)量最多，達(dá)到58個(gè)，分布在19個(gè)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)分析檢測(cè)到的關(guān)聯(lián)位點(diǎn)更加集中。

為了探討角果皮在硫苷合成和轉(zhuǎn)運(yùn)中的作用，對(duì)授粉后25天的角果皮硫苷含量也進(jìn)行SNP關(guān)聯(lián)分析（圖2-c）。結(jié)果表明，僅在A8、A9和C3染色體上各檢測(cè)到1個(gè)顯著關(guān)聯(lián)位點(diǎn)，遠(yuǎn)遠(yuǎn)少于種子的檢測(cè)結(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個(gè)高質(zhì)量GEM標(biāo)記，其中，25 834個(gè)位于A基因組，28 055個(gè)位于C基因組。以硫苷含量為因變量，GEM標(biāo)記為自變量，利用混合線性模型進(jìn)行關(guān)聯(lián)分析。以校正后的 Bonferroni值（P=0.05/53889=9.28×10-7）為顯著性標(biāo)準(zhǔn)，取負(fù)對(duì)數(shù)后的近似值為6.03。另外，在控制假陽(yáng)性率的情況下，利用多重假設(shè)檢驗(yàn)FDR=0.05來(lái)檢測(cè)可能由于Bonferroni值太過(guò)嚴(yán)格而遺漏的候選基因或關(guān)聯(lián)區(qū)間。

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

成熟種子硫苷含量共檢測(cè)到127個(gè)顯著GEM標(biāo)記，分布于除A4和A10染色體外的其他所有17條染色體上（圖3-b），其中在A2、A9、C2和C9染色體上出現(xiàn)了4個(gè)明顯的關(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個(gè)標(biāo)記，顯著性最高為-log10P=10.15；標(biāo)記最少的峰是GEMb-A2，僅有13個(gè)，且代表性標(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后角果皮硫苷含量共檢測(cè)到24個(gè)顯著性標(biāo)記，分布于 A1、A2、A8、A9、C2、C4和 C8染色體上（圖3-c）。C9染色體上有一個(gè)明顯的關(guān)聯(lián)峰，但未達(dá)Bonferroni顯著性閾值，因此采用假陽(yáng)性率FDR=0.05為新的閾值來(lái)篩選關(guān)聯(lián)區(qū)間，最終得到GEMc-A8、GEMc-A9、GEMc-C2和GEMc-C9四個(gè)關(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角果皮與成熟種子的硫苷含量共同檢測(cè)到 3個(gè)明顯重合的關(guān)聯(lián)區(qū)間，分別位于A9、C2和C9上（圖3和表5）。不僅在區(qū)間上存在重合，而且有22個(gè)顯著標(biāo)記在兩個(gè)性狀間也是相同的。而成熟種子與25 DAP種子的硫苷含量卻只檢測(cè)到1個(gè)共同的顯著標(biāo)記，25 DAP種子與25 DAP角果皮的硫苷則沒(méi)有檢測(cè)到共有的顯著標(biāo)記。

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

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

表6 SNP和GEM共同檢測(cè)到的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ù)測(cè)

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

3 討論

利用一個(gè)歐洲甘藍(lán)型冬油菜自然群體，首次開(kāi)展了不同組織器官硫苷含量動(dòng)態(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角果皮檢測(cè)到的顯著性GEM標(biāo)記大部分（92%）在成熟種子中能重復(fù)檢測(cè)到，暗示油菜中的硫苷也可能像擬南芥、蕪菁以及芥菜一樣存在由角果皮到種子的轉(zhuǎn)運(yùn)過(guò)程[10,34]。

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

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

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

4 結(jié)論

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

[1] 張永霞, 趙鋒, 張紅玲. 中國(guó)油菜產(chǎn)業(yè)發(fā)展現(xiàn)狀問(wèn)題及對(duì)策分析.世界農(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] 王漢中, 殷艷. 我國(guó)油料產(chǎn)業(yè)形勢(shì)分析與發(fā)展對(duì)策建議. 中國(guó)油料作物學(xué)報(bào), 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)因子及其對(duì)奶牛生產(chǎn)的影響. 中國(guó)草食動(dòng)物, 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)型油菜硫苷特性與對(duì)菌核病抗性關(guān)系. 湖南農(nóng)業(yè)大學(xué)學(xué)報(bào)(自然科學(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.