甘藍型油菜種子和角果皮中硫苷含量的動態(tài)變化及轉錄組關聯分析

田志濤 ，趙永國 ，Lenka Havlickova，He Zhesi，Andrea L Harper，Ian Bancroft，鄒錫玲，張學昆，陸光遠

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

0 引言

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

1 材料與方法

1.1 材料和樣品采集

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

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

1.2 硫苷檢測方法

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

1.3 轉錄組測序及SNP和GEM的檢測

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

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

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

續(xù)表1 Continued table 1

續(xù)表1 Continued table 1

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

1.4 SNP關聯分析

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

1.5 GEM關聯分析

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

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

2 結果

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

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

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

種子成熟后（約授粉后50 d），硫苷含量總體上進一步升高，平均為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種子的最高值。由此可見，硫苷的合成、運輸和積累是一個十分復雜的動態(tài)變化過程。

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

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

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

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

2.2 SNP關聯分析

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

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

圖2 硫苷含量的SNP關聯分析Fig. 2 SNP association mapping of glucosinolate contents

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

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

表4 成熟種子硫苷含量SNP關聯分析結果Table 4 SNP association peak for glucosinolate content in mature seeds

2.3 GEM關聯分析

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

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

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

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

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

圖3 硫苷含量的GEM關聯分析Fig. 3 GEM association mapping of glucosinolate contents

2.4 不同標記關聯分析比較

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

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

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

表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 候選基因的預測

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

表8 SNP和GEM關聯峰內篩選到的候選基因Table 8 Summary of candidate genes within association peaks with SNPs and GEMs

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

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

3 討論

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

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

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

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

4 結論

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

[1] 張永霞, 趙鋒, 張紅玲. 中國油菜產業(yè)發(fā)展現狀問題及對策分析.世界農業(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] 王漢中, 殷艷. 我國油料產業(yè)形勢分析與發(fā)展對策建議. 中國油料作物學報, 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] 申躍宇, 尹佩輝, 莫放. 菜籽餅(粕)的抗營養(yǎ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] 宋志榮, 官春云. 甘藍型油菜硫苷特性與對菌核病抗性關系. 湖南農業(yè)大學學報(自然科學版), 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.