小豆五出復葉pcl突變基因的精細定位和轉錄分析

劉 震 鄧 勇 趙超綺 星君宜 萬 平 楊 凱

（北京農學院，植物科學技術學院/農業應用新技術北京市重點實驗室，北京 102206）

小豆［Vigna angularis（Willd）Ohwi &Ohashi］，別名紅小豆、赤豆等，為一年生草本植物［1］，隸屬豆科（Leguminosae）豇豆屬（Vigna）［2］。小豆是我國主要的食用豆類之一，已有上千年的栽培歷史［3］，種質資源極為豐富，但對小豆農藝性狀的研究卻很薄弱。2015年，Yang 等［4］完成了小豆全基因組測序，從基因組層次加快了對小豆農藝性狀的研究分析。

葉片是小豆進行光合作用的重要器官，其在產量和品質的形成中發揮著重要作用［5］。多數豆科作物為復葉植物，其復葉形態能夠直接影響光合特性和生物量積累。前人研究表明，綠豆多小葉復葉種質增加了葉面積指數，能更有效地吸收光能，同時綠豆七小葉復葉品系促進了產量提高［6］；宗春美等［7-8］研究發現，大豆五小葉復葉品系牡5796-3具有較好的豐產性能。

突變體是開展作物遺傳研究和育種的基礎，具有優良性狀的突變體可以為育種、科研和生產提供優質的種質資源。宗春美等［8］創制了一個多小葉大豆品系牡5796-3，初步得到了多小葉性狀的遺傳規律。Chontira 等［9］也通過輻射獲得了綠豆多小葉性狀的突變體。Wang 等［10］發現了一種起源于野生大豆的多小葉突變體，可創制抗病材料［11］。這些研究豐富了豆科多小葉突變體的種質資源，為探究豆科植物復葉調控模式提供了優質材料。

目前在豌豆（Pisum sativum）突變體中發現了uni單葉突變體［12］、cri異位小葉突變體［13］以及能夠降低葉片復雜程度的stp突變體［14］。在蒺藜苜蓿（Medicago truncatula）中發現了sgl1單葉突變體［15］、palm1五出復葉突變體［16］、pinna1五出復葉突變體［17］。在綠豆（Vigna radiata）中發現了un單葉突變體［18］。在百脈根（Lotus japonicus）中發現了能夠導致5 片小葉有不同程度融合，嚴重時5片小葉融合為1片葉的uml突變體［19］。但目前與小豆復葉發育相關的突變體研究還鮮有報道。

因此，北京農學院小豆課題組（本課題組）前期通過使用電子束射線在北京輻射中心對京農6 號供試種子進行輻射，在其輻射后代中篩選出一株五小葉復葉個體。該突變體較常規三小葉復葉小豆品種多了2 個小葉，經多代自交，獲得了可穩定遺傳的五出復葉葉型材料。本研究以小豆五出復葉突變體為材料，對小豆復葉發育調控基因進行遺傳分析與圖位克隆，確定五出復葉發育調控候選基因，并對野生型和突變體小豆進行轉錄組測序，分析突變性狀關鍵候選基因可能的調控途徑與功能，以期對小豆復葉發育的遺傳調控提供研究方向。

1 材料與方法

1.1 試驗材料

三出復葉小豆栽培品種GM635 與五出復葉突變體EB601 由北京農學院小豆課題組提供：使用300 Gy計量電子束處理小豆栽培品種京農6號，在后代中發現一株五出復葉突變體（圖1），經過5 代自交獲得遺傳穩定的五出復葉材料（EB601），復葉性狀基因暫定名為pcl（pentafoliate compound leaf）。

圖1 三出復葉和五出復葉的表型特征Fig.1 Phenotypic characteristics of ternately and pentafoliate leaf

構建EB601×GM635 雜交組合，F1均為三出復葉表型，F2出現性狀分離，收獲F2單株，選用F2∶3株系用于遺傳分析和基因初定位，F3∶4分離群體用于基因精細定位。

1.2 農藝性狀分析

小豆農藝性狀調查按照《小豆種質資源描述規范和數據標準》［20］進行。

1.3 突變體遺傳分析與基因定位

1.3.1 突變體遺傳分析 構建EB601×GM635雜交組合，調查F2∶3株系農藝性狀并記錄三出復葉和五出復葉表型，進行突變體遺傳分析。

1.3.2 小豆基因組DNA 提取 采用十六烷基三甲基溴化銨（cetyltrimethylammonium bromide，CTAB）法［21］提取小豆基因組DNA。

1.3.3pcl的初步定位 利用自主開發的小豆1 572對簡單重復序列（simple sequence repeats，SSR）標記，對親本和F2∶3株系差異性狀個體基因組DNA 分別混池，并進行PCR 擴增篩選。找出與目標基因連鎖的標記，再通過F2∶3家系對標記進行群體驗證。PCR 反應體系共15 μL：5 U·μL-1TaKaRa Taq 0.25 μL，10 μmol·L-1正反引物各1 μL，1 mmol·L-1dNTPs 1.5 μL，1.5 mmol·L-1Mg2+1 μL，50 ng 小豆基因組DNA 2 μL，10×PCR Buffer 1.5 μL，ddH2O 6.75 μL；PCR 擴增程序：95 ℃預變性5 min，95 ℃變性30 s，58 ℃退火30 s，72 ℃延伸20 s，72 ℃終延伸10 min，35個循環。

1.3.4pcl的精細定位 提取F3∶4家系中個體和親本的基因組DNA；三出復葉親本和五出復葉親本各取5 株DNA 樣本分別等量混池，F3∶4家系中三出復葉個體和五出復葉個體各取50 株DNA 樣本分別等量混池。通過高通量測序技術（“next-generation” sequencing technology，NGS）和集群分離分析法（super bulked segregant analysis，Super BSA）以及特異位點擴增測序技術（specific-locus amplified fragment sequencing，SLAF-seq）對親本混池和差異個體混池進行測序分析，開發差異位點?；跀祿Y果使用Oligo 7.0 軟件開發SSR 與單核苷酸多態性（single nucleotide polymorphisms，SNP）標記，構建三出復葉個體基因組DNA混池（10株）和五出復葉個體基因組DNA混池（10株），對親本和2個差異性狀個體基因組DNA混池進行引物篩選。反應體系與初定位相同，使用瓊脂糖凝膠電泳檢測PCR產物。通過SSR和SNP標記對初定位區間進行加密，縮小定位區間。

1.4 染色體步行

參考小豆基因組序列［4］，對定位區間進行染色體步行。使用Oligo 7.0 軟件設計克隆引物，引物合成由生工生物工程（上海）股份有限公司完成。使用2×Phanta?Max Master Mix（諾唯贊，南京）進行PCR 擴增，瓊脂糖凝膠電泳回收PCR 產物，由深圳華大基因科技有限公司完成測序，使用啟動區域元件分析網（http：//bioinforma-tics.psb.ugent.be/webtools/plantcare/html）分析克隆結果中啟動區域產生的突變位點。

1.5 轉錄組分析

1.5.1 RNA 提取和質量檢測 采用TRIzol 試劑盒（賽默飛，美國）提取小豆三出復葉和小豆五出復葉幼葉樣品的總RNA。通過Agilent 2100 Biognlyzer生物分析儀（安捷倫，美國）和NanoDrop 紫外分光光度計（賽默飛，美國）檢測提取RNA 樣品的質量、純度及完整性，OD260/280在2.0～2.2 之間，RNA 分子完整數在6.6～8.5之間，可用于建庫和測序。

1.5.2 cDNA文庫構建、測序 選取3個生物學重復材料構建RNA-seq 文庫。質量檢測合格后，用Oligo（dT）磁珠富集樣品中mRNA。加入打斷試劑（fragmentation buffer）將mRNA 打斷成短片段，反轉錄合成雙鏈cDNA。利用AMPure XP beads 磁珠純化雙鏈cDNA，末端修復加A尾并連接測序接頭。最后PCR擴增構建cDNA 文庫。用Agilent 2100 Biognlyzer 生物分析儀和ABI StepOnePlus RealTime PCR System（賽默飛，美國）對文庫的插入片段范圍進行檢測，對文庫濃度進行定量，質檢合格后，使用Illumina HiSeq 測序儀（Illumina，美國）進行測序。

1.5.3 原始數據組裝 對測序所得原始數據進行過濾處理，去除包含測序接頭的reads；去除未知堿基N含量大于10%的reads；去除低質量堿基（Q＜5）含量大于50%的reads。統計clean reads 的數量、總長度、Q20、Q30、GC%等。將過濾后數據與Yang 等［4］的小豆參考基因組進行比對。

1.5.4 差異表達基因的篩選和功能注釋 采用基于負二項分布模型的算法（DEGseq）進行差異表達分析。對差異基因進行篩選的閾值一般為Log2（Fold Change）＞1 且P＜0.05。通過京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析，對差異表達基因進行通路分析，統計每個通路中包含的差異表達基因的個數，通過超幾何檢驗找出顯著富集的通路。

1.5.5 實時熒光定量PCR（quantitative real-time PCR，qRT-PCR）驗證 采用TB Green?Premix Ex Taq? II試劑盒（TaKaRa，日本），采用CFX96 實時熒光定量PCR 系統（Bio-Rad，美國）進行qRT-PCR 分析，qRTPCR 總體積為20 μL：TB Green Premix Ex Taq II（2X）9 μL，10 μmol·L-1正反引物各1 μL，100 ng cDNA 2 μL，ddH2O 7μL。反應條件：95 ℃預變性30 s；95 ℃變性10 s，58 ℃退火15 s，39 個循環；溶解曲線72 ℃ 30 s。以小豆β-Actin基因（Va03G001640）作為內參基因，每個樣品分別設3次生物學重復和3次技術重復，以2-ΔΔCT法計算基因的相對表達量，并計算三出復葉野生型和五出復葉突變體之間的差異倍數Log2（Fold Change）。

1.6 數據分析

1.6.1 農藝性狀分析 使用SPSS 19.0 和Excel 2007軟件，對調查數據進行統計分析。

1.6.2 定位區間連鎖圖 使用Mapmaker 3.0 軟件［22］構建基因連鎖圖，設置LOD=3.0，最大圖距50.00 cM。

2 結果與分析

2.1 農藝性狀分析

對EB601×GM635組合F2：3株系個體的株高、莖粗、主莖節數、單株莢數、莢長、莢寬、單莢粒數、單株產量、百粒重共9 個主要農藝性狀進行調查統計分析，結果如表1所示。與三出復葉個體相比，五出復葉個體株高變高，莖變粗，主莖節數增多，豆莢和籽粒變大，單株莢數和單莢粒數增多，單株產量和百粒重提高；其中株高、莖粗、莢寬、單莢粒數和百粒重呈顯著性差異（P＜0.05），單株產量呈極顯著性差異（P＜0.01）。

表1 EB601×GM635 F2∶3株系農藝性狀統計Table 1 Statistics of interesting agronomic traits of EB601×GM635 F2∶3 segregating population

對EB601×GM635 組合F2：3株系間9 個主要農藝性狀進行了相關性分析。結果表明，在三出復葉和五出復葉表型之間，株高、莖粗、主莖節數、單株莢數之間的相關性顯現出較好的一致性，均為極顯著相關；而株高與莢長、單莢粒數、百粒重之間，單株莢數與莢長之間，單株產量與莢寬、百粒重之間，均未表現出共同的一致性。單株產量與各農藝性狀之間，只有五出復葉個體在莢寬上表現為極顯著相關（表2）。

表2 EB601×GM635 F2∶3各農藝性狀相關性分析Table 2 Correlation analysis of agronomic traits in EB601×GM635 F2∶3 segregating population

綜上，五出復葉突變基因對農藝性狀有較大的影響，且會導致農藝性狀之間的相關性發生變化。

2.2 小豆五出復葉調控基因的遺傳分析

2016年種植的EB601×GM635 雜交組合F2群體單株收獲，調查F2中三出復葉和五出復葉表型情況。2017年共種植了12 個EB601×GM635 雜交組合F2三出復葉個體的F2：3株系和3個F2五出復葉個體的F2∶3株系。3 個F2五出復葉個體的F3株系中個體沒有復葉性狀分離，均為五出復葉；12個F2三出復葉個體的F3株系中有7個株系出現復葉性狀分離，對分離株系群體中個體復葉性狀分離比（三出復葉個體數：五出復葉個體數）進行卡方檢驗，均符合3∶1的假設（表3）。

表3 EB601×GM635雜交組合F2∶3株系個體復葉性狀分離統計Table 3 Segregation of pentafoliate and trifoliolate among F2∶3 family of EB601×GM635 populatuon

F2五出復葉個體的F2∶3株系中個體沒有復葉性狀分離，F2三出復葉個體的F2∶3株系中個體出現復葉性狀分離現象，且三出復葉個體和五出復葉個體的分離比符合3∶1，表明小豆五出復葉個體由一對隱性基因控制，且符合孟德爾遺傳。

2.3 五出復葉調控候選基因定位

隨機選取GM635×EB601 F2：3株系三出復葉和五出復葉個體各10 株，將DNA 等量混合為三出復葉DNA混池和五出復葉DNA 混池。利用本課題組自主開發的1 572 對小豆基因組SSR 標記，篩選親本和F2∶3差異個體DNA 混池間的差異標記，僅有BAgs369、BAgs802、BAgs1061這3對SSR標記擴增產物電泳帶型滿足pcl基因標記特征。使用這3 個SSR 標記進行群體驗證，確定這3 個SSR 標記擴增位點與pcl基因存在連鎖關系，標記均位于小豆參考基因組第5 號染色體的末端。根據參考基因組序列在定位區間內的SSR位點，開發了86 個SSR 標記。pcl基因定位在標記C5S207-5 和標記BAgs1061 之間13.4 cM 的候選區域內（表4、圖2），標記C5S207-5 和BAgs1061 分別距pcl基因距離9.6和3.8 cM。

使用NGS 和Super-BSA 對GM635×EB601 F3∶4株系的三出復葉和五出復葉DNA 混池和親本混池進行測序分析。獲得的原始數據產量（Raw reads）平均值為15 500 821，過濾后有效數據量（Clean base）平均值為4.64 Gb，Clean reads 平均值為15 454 177，有效堿基率（Effective rate）平均值為99.70%，樣品Q20 平均值為96.93%，GC 含量占比平均值為35.70%（表5）。上述結果表明，本研究測序質量符合要求?；赟LAF-seq技術發現，連鎖區域內共有132 個SNP 位點，等距開發設計了23對SNP引物，并對親本和F3：4群體254個單株進行驗證，最終篩選出3 個SNP 標記88601FWR、15106FRM 和89739FMR 與pcl基因連鎖（表4、圖2）。進而將調控五小葉復葉的基因pcl定位到SNP 標記89739FMR 和SNP 標記15106FRM 之間，標記88601FWR、15106FRM 和89739FMR 距pcl基因距離分別為1.0、1.2和0.4 cM（圖3）。

表5 Super-BSA 重測序數據統計Table 5 Basic date for samples by super-BSA

圖2 部分EB601×GM635家系分離群體個體開發引物擴增電泳圖譜Fig.2 Electrophoretic map of developed primer amplification for EB601×GM635 segregating populations individuals

表4 開發的部分標記和復葉發育相關基因表達分析引物Table 4 Some developed markers and quantitative primers of compound leaf development related genes

2.4 五出復葉定位區間染色體步行

依據京農6號參考基因組，對基因定位區間1.6 cM內的基因進行分析。連鎖標記的遺傳與物理位置排列順序一致，對應小豆參考基因組5號染色體上約500 kb的區間，共有27個基因（表6）。通過染色體步行發現，定位區間內距VaTMK3基因前960 bp 位置的啟動區存在一個突變堿基（A 變C）位于TATA-box 元件上，該突變位點分析與性狀共分離（圖3）。VaTMK3編碼TMK類蛋白，參與細胞擴張、生長素轉運。葉片細胞的增殖擴張與TMKs密切相關，TMKs通過調節生長素來控制植物的生長發育。

表6 500 kb定位區間基因注釋Table 6 Gene annotation in 500 kb localized interval

2.5 轉錄組分析

2.5.1 基因表達量分析 三出復葉和五出復葉幼葉樣本中檢測到基因數量在24 800～25 529個之間，對樣本中基因表達量進行統計，結果顯示每個樣本中有2 116（8.44%）～4 552（18.63%）個基因的FPKM（Fragments Per Kilobase per Million）＞30，8 999（36.89%）～10 166（40.47%）個基因的FPKM 在5～30之間，占比最高。三出復葉和五出復葉在幼葉中有5 260個差異表達基因，相對于三出復葉幼葉，五出復葉幼葉中有3 411個基因表達上調，有1 849個基因表達下調（圖4）。

圖4 差異表達基因分析Fig.4 Analysis of differentially expressed genes

2.5.2 差異表達基因Pathway 富集分析 KEGG 富集結果顯示，在幼葉中差異表達基因被富集到133條通路中，其中顯著富集的有14 條（P＜0.05），主要包括植物MAPK 信號通路（MAPK signaling pathway-plant，ko04016）、植物病原互作（plant-pathogen interaction，ko04626）（表7）。植物MAPK 信號轉導通路參與植物的生長發育，與植物激素的生物合成、運輸和信號轉導之間存在緊密聯系。說明小豆五出復葉調控模式與植物激素合成運輸信號轉導有關。

表7 差異表達基因KEGG富集Table 7 KEGG enrichment of differentially expressed genes

2.5.3 五出復葉中相關生長素調控基因差異表達轉錄組測序與qRT-PCR 數據分析結果發現，相對于三出復葉中的表達情況，生長素信號傳導相關基因在五出復葉中顯著或極顯著差異表達，如VaTMK3在幼葉中表達極顯著上調；參與生長素轉運與介導生長素梯度形成的VaPIN1a（auxin efflux carrier component 1-like）和VaPIN1b（auxin efflux carrier component 1-like）基因在五出復葉幼葉中均極顯著上調表達；同時在五出復葉幼葉中參與質子驅動的生長素流入的載體蛋白編碼基因VaLAX1（auxin transporter-like protein 1）、VaLAX3（auxin transporter-like protein 3）和VaLAX5（auxin transporter-like protein 5）均上調表達，這類基因可介導從發育中的葉片（生長素生物合成部位）到尖端的生長素梯度的形成；編碼生長素轉運抑制蛋白（transport inhibitor response 1-like protein）的基因Va07G049160和Va03G081560也呈現上調表達（圖5）。說明生長素濃度與梯度分布可能參與小豆五出復葉突變體發育調控。

圖5 轉錄組測序與qRT-PCR分析基因相對表達Fig.5 The relative expression of key genes by transcriptome sequencing and qRT-PCR analysis

2.5.4 五出復葉中葉片極性調控相關基因差異表達 轉錄組測序與qRT-PCR 數據分析發現，在五出復葉幼葉中能夠調節葉片極性的YABBY 基因家族成員VaYAB5（axial regulator YABBY5）、VaYAB1（axial regulator YABBY1）和VaYAB4（axial regulator YABBY4）基因表達量顯著上調（圖5）。說明小豆YABBY基因家族基因過表達可能參與五出復葉突變體葉片極性調控。

3 討論

本研究通過小豆五出復葉性狀遺傳分析發現，三出復葉性狀為顯性性狀。五出復葉小豆較三出復葉小豆生長速度快，植株高，生長茂盛，說明pcl基因對植株的生長影響程度較大。在小豆五出復葉pcl的作用下，株高變高，莖變粗，而且主莖節數增多，說明小豆五出復葉突變體對小豆生長發育具有積極作用。由上述結果可知，研究小豆五出復葉突變體對闡述小豆復葉發育機理和小豆葉型遺傳改良具有重要意義。

前人研究發現TMKs 基因家族通過調節細胞擴張和細胞增殖來協調植物生長［23］。同時TMKs還參與生長素信號轉導調控途徑TMK1-IAA32/34-ARFs 和TMK1/ABP1-ROP2/6-PINs。擬南芥中，TMKs的功能與在葉片起始和發育的早期增殖階段直接調控細胞的生長有關［23］；TMKs 家族基因表達下降導致擬南芥對生長素的敏感性下降，所以TMKs 家族似乎通過調節細胞擴張和增殖來參與植物葉片的生長，并作為生長素信號傳導的一個組成部分影響植物中的生長素轉運［23］。本研究對小豆復葉發育調控基因進行遺傳分析與圖位克隆，確定五出復葉發育調控候選基因為VaTMK3。VaTMK3屬于TMKs 家族，編碼的TMK 蛋白與生長素信號傳導以及細胞生長發育密切相關［23］。故推測五出復葉突變體VaTMK3表達的上調可能會影響生長素信號轉導。進一步分析發現，五出復葉突變體幼葉中參與質子驅動的生長素流入的載體蛋白基因VaLAX1、VaLAX3、VaLAX5表達均上調；同時編碼生長素轉運抑制蛋白的基因Va07G049160和Va03G081560基因也都上調表達，說明VaTMK3基因表達量的改變可能對調控生長素轉運基因產生影響。

前人對番茄、碎米芥的復葉研究發現，復葉原基的葉邊緣分裂區域中，存在高濃度生長素，其對小葉生長、葉原基發生均產生了影響，且單葉和復葉里生長素濃度變化與PIN（PIN-FORMED）基因密切相關［24］。已報道發現百脈根中PIN1同源基因UML的突變體會導致小葉原基融合，頂葉出現“杯子”、“雨傘”形態［19］。TMKs能夠通過調控參與介導植物生長素運輸的基因PIN（PIN-FORMED）來影響生長素的分布［25］。本研究在五出復葉突變體小豆幼葉中發現VaPIN1a、VaPIN1b均上調表達，表明五出復葉中VaTMK3表達量的上調可能調控了VaPIN1a和VaPIN1b的表達，從而影響了小豆復葉形態。

當前研究還發現，生長素外排轉運體PIN3 和生長素內流載體LAX3 的連續作用能夠間接激活MAPK 信號通路中MPK3/MPK6 級聯反應，從而促進CWR基因的表達，導致側根的產生［26］。而小豆五出復葉中MAPK 信號通路差異顯著，推測小豆額外小葉的產生可能與側根產生的模式相似，生長素的積聚可能起到關鍵作用。

在五出復葉突變體幼葉中，YABBY 家族基因VaYAB1、VaYAB4和VaYAB5表達量明顯上調。有報道發現YAB1能夠參與器官發生過程中遠軸細胞形態的建立并調節胚胎芽頂端分生組織（shoot apical meristem，SAM）發育的起始［27］，且轉錄因子YAB1在擬南芥中影響葉發育；YAB4對胚珠外珠被的形成和背面-正面不對稱生長至關重要，同時促進近軸細胞特性［28］；YAB5調節胚芽頂端分生組織的起始發育［29］。在葉器官發育起始時YAB1和YAB5在葉片遠軸面表達，促進葉極性的確立［30］。在大豆（Glycine max）中，YAB1類基因GmFILa在擬南芥中過表達引起葉片近軸端表現為部分遠軸化，導致葉片卷曲，生育期延長，同時影響自由態生長素含量［31］。由于生長素活性在葉片的近遠軸間的差異對葉片形態發育產生影響［32］，推測VaYAB1、VaYAB4和VaYAB5基因表達可能參與生長素濃度梯度分布的變化，進而影響小豆復葉發育過程，但YABBY 家族基因是否直接參與小豆復葉發育還需進一步驗證。

綜上所述，小豆復葉發育調控可能是由候選基因VaTMK3通過影響生長素跨膜轉運相關基因VaPIN1a、VaPIN1b、VaLAX1、VaLAX3、VaLAX5，調控小豆葉片中生長素濃度梯度，導致分裂帶發生生長素積累，從而調控異位小葉起始發育；同時，生長素濃度改變會導致YABBY 基因家族中VaYAB1、VaYAB4 和VaYAB5表達上調，致使小豆葉片遠軸極性發生變化，即在YaTMK3和YABBY家族基因共同作用下導致五出復葉產生。

4 結論

本研究通過輻射誘變獲得小豆五出復葉突變體pcl，相比于野生型小豆，其農藝性狀均有所提高。遺傳分析結果表明，五出復葉突變體pcl受一對隱性基因控制，定位區間內確定候選基因為VaTMK3，該基因參與生長素信號轉導，其在突變體pcl幼葉中顯著上調表達。突變體pcl幼葉中參與生長素跨膜轉運載體蛋白基因VaPIN1a、VaPIN1b、VaLAX1、VaLAX3、VaLAX5均上調表達。由此推測VaTMK3可能是通過影響生長素濃度梯度導致小豆復葉的產生改變，以及YABBY家族可能參與了小豆復葉的發育過程。