孫豐磊,任姣姣,雷 斌,高文偉,曲延英
(1.新疆農(nóng)業(yè)大學(xué)農(nóng)學(xué)院,烏魯木齊 830052;2 新疆農(nóng)業(yè)科學(xué)院核技術(shù)生物研究所,烏魯木齊 830052)
【研究意義】不同棉花品種的不同生長發(fā)育階段對干旱脅迫響應(yīng)不同。棉花花鈴期若受到干旱脅迫會影響棉花產(chǎn)量。以現(xiàn)有棉花材料為基礎(chǔ),篩選抗旱性強(qiáng)材料,挖掘、研究棉花抗旱的相關(guān)基因,對提高棉花的抗旱性及加快培育抗旱棉花新品種的進(jìn)程有重要意義?!厩叭搜芯窟M(jìn)展】作物抗旱性是一個復(fù)雜的數(shù)量性狀,傳統(tǒng)的定位方法耗時長、成本高[1]。混池測序是一種基于二代測序以集群分離方法為主的技術(shù),應(yīng)用于基因定位和克隆等[2]。BSA(集群分離方法)應(yīng)用在萵苣的抗霜霉病研究中,篩選到3個RAPD分子標(biāo)記[3],3個分子標(biāo)記與抗性基因Dm5/8緊密連鎖[4]。通過選擇具有極端性狀或具有明顯表型差異的親本來構(gòu)建群體,從子代群體中根據(jù)目標(biāo)性狀的特點選擇具有極端差異子代個體,構(gòu)成兩個極端性狀的子代DNA混池,通過測序分析獲得與性狀相關(guān)的分子標(biāo)記,定位相關(guān)基因。由于作為兩個親本是在目標(biāo)性狀中選擇明顯差異的個體,因此,在子代的兩個極端混池中,只有在目標(biāo)性狀基因的區(qū)域中有差異,其它區(qū)域的差異較小或基本沒有差異,增加了定位的準(zhǔn)確性和可靠性。混池測序(BSA-seq)早期應(yīng)用在基因組較小的物種中,并且目標(biāo)基因在該物種的突變體中。Schneeberger等[5]利用混池測序成功定位到導(dǎo)致擬南芥葉片淺綠色突變體的基因。在擬南芥中還定位到配子突變導(dǎo)致不育的基因[6]以及突變導(dǎo)致開花延遲的基因[7]。而此方法在水稻中也得到應(yīng)用,Abe等[8]在EMS誘變的突變體水稻中將白化和半矮稈定位在某一區(qū)間內(nèi),進(jìn)而對目標(biāo)進(jìn)行了定位和克隆。Takagi等[9]在EMS誘變的水稻突變體中,定位到一個耐鹽堿基因,通過MutMap 的方法克隆該基因并培育出耐鹽堿的品種。陳忠明等[10]結(jié)合BSA分析法和RFLP的方法,在可育與半不育的極端混池中發(fā)現(xiàn)了2個標(biāo)記(RG213、G200)與親和基因連鎖。陸光遠(yuǎn)等[11]利用油菜的回交群體篩選到兩個AFLP標(biāo)記,2個標(biāo)記與細(xì)胞核雄性不育基因緊密連鎖。李傳友等[12]結(jié)合BSA和AFLP分析法篩選到光敏核不育基因相關(guān)的多態(tài)性產(chǎn)物。張顯等[13]在西瓜隱形和雄性不育和可育的植株中,通過RAPD-BSA方法篩選出s1167的引物可以區(qū)分該兩種株系。Wang等[14]通過BSA分析的方法在水稻中篩選到9個分子標(biāo)記,利用9個分子標(biāo)記最終定位到tms5基因。Geng等[15]利用BSA分析法對油菜千粒重基因進(jìn)行定位,篩選到4個與千粒重相關(guān)的候選基因。Hong等[16]通過該方法篩選到與甜味相關(guān)的候選基因1個,酸味相關(guān)的候選基因8個。Zheng等[17]通過混池測序,共篩選到4個抗稻瘟病候選基因。而Zhang等[18]通過混池關(guān)聯(lián)分析的方法研究黃瓜抗白粉病相關(guān)基因,共篩選到6個相關(guān)的候選基因。Chen等[19]和Feng等[20]通過BSA測序分析,再結(jié)合SSR標(biāo)記,定位到與抗稻瘟病相關(guān)的基因。Qin等[21]在大麥白化性狀方面利用集群分離分析的方法,結(jié)合SSR分子標(biāo)記,定位6個有關(guān)該性狀的候選基因。杜威世等[22]在海島棉品種和陸地棉品種的雜交F2代中,結(jié)合SSR分子標(biāo)記和BSA法定位到3個與抗棉花黃萎病相關(guān)的基因標(biāo)記。Hayes等[23]研究結(jié)果表明,在大豆中通過BSA法篩選到4個Rsv1相關(guān)的分子標(biāo)記,該4個分子標(biāo)記與抗花葉病相關(guān)?!颈狙芯壳腥朦c】棉花抗旱性是一個復(fù)雜的數(shù)量性狀,傳統(tǒng)的定位方法耗時長、成本高?;斐販y序是一種基于二代測序以集群分離方法為主的技術(shù),可以有效用于棉花復(fù)雜性狀的遺傳定位。【擬解決的關(guān)鍵問題】篩選出抗旱性強(qiáng)和抗旱性敏感的材料,構(gòu)建抗旱性強(qiáng)和抗旱性弱的兩個極端樣本混池,采用BSA混池測序技術(shù)測序分析親本及子代混池,定位篩選抗旱基因位點候選區(qū)域,為挖掘棉花相關(guān)的抗旱基因奠定基礎(chǔ)。
選用由石遠(yuǎn)321為父本和奎85-174為母本構(gòu)建的150份材料(重組自交系子代群體2年2點田間抗旱表型鑒定),篩選得到抗旱性強(qiáng)的20份極端材料和20份水分敏感的極端材料,分別提取極端材料的葉片DNA,檢測濃度以及純度,等量混合,分別構(gòu)成極端抗旱池與極端敏旱池,提取2個親本的葉片DNA。
樣本混池建庫和測序由北京百邁客生物科技有限公司完成。按照lllumina公司執(zhí)行標(biāo)準(zhǔn),首先檢測DNA樣品后,進(jìn)行超聲破碎[24],修復(fù)、拼接、擴(kuò)增,并構(gòu)建完成測序所需要的文庫。經(jīng)過質(zhì)檢合格后的文庫用測序分析(參考陸地棉基因組(版本號為v2.1))。
1.2.1 SNP變異位點的檢測
SNP變異位點檢測是通過GATK[25-26]軟件工具包完成。根據(jù)Clean Reads在參考基因組的定位結(jié)果,通過去重復(fù)(Mark Duplicates)以及GATK進(jìn)行局部重比對(Local Realignment)、堿基質(zhì)量值校正(Base Recalibration)等處理,再使用GATK檢測單核苷酸多態(tài)性(Single Nucleotide Polymorphism,SNP),得到SNP位點集。
根據(jù)Clean Reads在參考基因組上的定位結(jié)果,小片段的插入與缺失有基因組之間的檢測確定。但同樣反映了其差異的Small InDel變異的數(shù)量較少,其中InDel時會引起移碼突變[27],如果移碼突變區(qū)在編碼區(qū)將導(dǎo)致基因功能上的改變。
1.2.2 關(guān)聯(lián)分析及候選區(qū)域的基因注釋
對SNP過濾后,進(jìn)行關(guān)聯(lián)分析。在SNP關(guān)聯(lián)分析中主要采用歐式距離(Euclidean Distance,ED)和SNP-inde 2種分析方法,將2種分析方法共同關(guān)聯(lián)到的區(qū)間作為抗旱相關(guān)的主要區(qū)間。
其中歐氏距離分析方法(Euclidean Distance,ED)是通過分析測序數(shù)據(jù)在混池間確定差異顯著的標(biāo)記,確定與目標(biāo)性狀關(guān)聯(lián)的區(qū)域。在目標(biāo)性狀相關(guān)位點上除了2個子代混池與親本間存在顯著差異外,其它位點不存在差異或趨于基本一致,通過分析混池間差異的SNP位點,計算ED值,在實驗中取ED的2次方進(jìn)行擬合。
1.2.3 熒光定量驗證候選基因的差異表達(dá)
將10份抗旱性極端材料在實驗室內(nèi)進(jìn)行土培(土壤為黑土∶蛭石為2∶1混合)培養(yǎng)至三葉期時,進(jìn)行土壤水分脅迫試驗,用于熒光定量測定基因表達(dá)量。
統(tǒng)計混池間基因型頻率差異的分析方法應(yīng)用在標(biāo)記關(guān)聯(lián)分析中稱為SNP-index算法[28]。是通過Δ(SNP-index)統(tǒng)計分析混池間基因型頻率間的顯著差異。選擇性狀相關(guān)的區(qū)域是在標(biāo)記的染色體上通過ΔSNP-index值進(jìn)行擬合分析后關(guān)聯(lián)閾值以上的區(qū)域。
通過BLAST軟件對候選基因在多個數(shù)據(jù)中(GO[20]、KEGG[21]、NR[22])比對詳細(xì)注釋。
研究表明,測序得到原始reads為119 060 015、121 194 054、190 811 397、223 420 448,得到clean reads與原始reads基本一致,總數(shù)據(jù)量為197.65 Gbp,其中Q30(%)在93%以上,并且GC含量在35.43%~35.73%,共得到142.23 Gbp的Clean reads,平均每個樣品測序深度15.47X,最終的平均比對效率為99.045%,達(dá)到18.75X平均覆蓋深度和99.03%的基因組覆蓋度。表1,表2
表1 樣品測序數(shù)據(jù)質(zhì)量統(tǒng)計
表2 與參考基因比對
研究表明,親本之間共獲得1 145 784個SNP,其中非同義突變的SNP共13 753個;混池之間共獲得588 221個SNP,其中引起非同義突變的SNP共6 186個。親本之間共獲得221 481個Small InDel;混池之間共獲得134 957個Small InDel。圖1
圖1 樣品間SNP(左)和InDel(右)Venn圖
2.3.1 SNP檢測
研究表明,共得到2 561 910個SNP位點,再次過濾后最終得到高質(zhì)量SNP位點共540 596個。
取所有位點的擬合值作為關(guān)聯(lián)閾值,得到ED值為0.17,共得到10個區(qū)域,總長度為39.30 Mb,共包含979個基因,其中包含非同義突變SNP位點的基因共112個。圖2,表3
注:橫坐標(biāo)為染色體名稱,彩色的點代表每個SNP位點的ED值,黑色的線為擬合后的ED值,紅色的虛線代表顯著性關(guān)聯(lián)閾值,下同
表3 ED算法關(guān)聯(lián)候選區(qū)域
擬合后的ΔSNP-index值為0.31,共得到5個區(qū)域,總長度為6.12 Mb,共包含86個基因,在所有基因非同義突變的共有27個。圖3,表4
圖3 ΔSNP-index關(guān)聯(lián)值在染色體上分布
表4 SNP-index算法關(guān)聯(lián)候選區(qū)域
2.3.2 候選區(qū)域篩選及候選區(qū)域基因注釋
研究表明,共篩選到3個候選區(qū)域,總長度為6.12 Mb,共包含86個基因。表5
表5 ED和SNP-index兩種算法關(guān)聯(lián)候選區(qū)域
親本間存在非同義突變的SNP共43個,混池間存在非同義突變的SNP共20個,在候選區(qū)域內(nèi)共注釋到83個基因,而注釋到的非同義突變基因在親本間共有27個。其中在GO數(shù)據(jù)庫中注釋到49個基因,而在KEGG數(shù)據(jù)庫中注釋到23個基因。
候選區(qū)域內(nèi)每個組分中基因數(shù)量不同,其中大多數(shù)被劃分到生物過程(Biological process)中。前5的功能分類分別是:在細(xì)胞組分中為質(zhì)體部分、葉綠體部分、細(xì)胞器膜、轉(zhuǎn)移酶復(fù)合體和葉綠體;在分子功能中為RNA導(dǎo)向的DNA聚合酶活性、核糖二磷酸羧化酶活性、蛋白質(zhì)結(jié)構(gòu)域特異性結(jié)合、雙加氧酶活性、碳水化合物磷酸酶活性;生物過程中主要為DNA整合、依賴RNA的DNA復(fù)制、糖醇代謝過程、細(xì)胞壁組織或生物合成、有機(jī)羥基化合物代謝過程。圖4,表6
圖4 候選區(qū)域內(nèi)基因GO注釋聚類
表6 候選區(qū)域內(nèi)基因GO注釋
得到的代謝通路中基因注釋分布最多的在代謝作用和遺傳信息處理中,代謝作用通路中注釋基因最多的是內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工和RNA聚合酶;在遺傳信息處理的代謝通路上苯丙烷生物合成和氨基酸生物合成中注釋到的基因種類最多。圖5
2.3.3 候選基因的篩選
研究表明,參與細(xì)胞組成的相關(guān)基因與棉花的抗旱性相關(guān),包括細(xì)胞核(GO:0005634),葉綠體(GO:0009507),過氧化物酶體(GO:0005777),RNA聚合酶Ⅲ復(fù)合物(GO:0005666)等,GH_A07G0705、GH_A07G0701、GH_D07G1930、GH_D07G1932、GH_D07G1934、GH_A07G0708、GH_A07G0695、GH_A07G0703等8個基因參與調(diào)控棉花的抗旱性。
鋅離子結(jié)合(GO:0008270)、RNA定向的DNA聚合酶活性(GO:0003964)、蛋白質(zhì)結(jié)構(gòu)域特異性結(jié)合(GO:0019904)、甲基轉(zhuǎn)移酶活性(GO:0008168)等參與到干旱脅迫相關(guān)反應(yīng)中。GH_A07G0695、GH_A10G1602、GH_D07G1927、GH_A10G1623、GH_A10G1601等24個基因與干旱脅迫調(diào)控相關(guān)。
一些參與干旱脅迫調(diào)控的過程,如代謝過程(GO:0008152),響應(yīng)脫落酸反應(yīng)過程(GO:0009737),有機(jī)物的運輸過程(GO:0071702),還有氧化還原過程(GO:0055114)等均參與干旱脅迫的生物調(diào)控過程中。預(yù)測GH_A07G0696、GH_A07G0700、GH_A07G0695、GH_A07G0703、GH_A07G0695、GH_A07G0696、GH_A07G0705等24個基因與干旱脅迫相關(guān)。表7
篩選出與前期5個棉花抗旱性評價關(guān)鍵性狀指標(biāo)相關(guān)的15個候選基因,分布在3個候選區(qū)域內(nèi)。表8
表8 抗旱候選基因及注釋
研究表明,在水分脅迫達(dá)到輕度(土壤含水下降了20%)、中度(土壤含水下降了40%)、重度(土壤含水下降了50%)時,在脅迫程度增加時,GhalkB隨著脅迫程度的增加都有一個上調(diào)表達(dá)的趨勢,在石遠(yuǎn)321中該基因在土壤含水量下降了10% 時表達(dá)量最高,隨著脅迫程度的增加,該基因的表達(dá)量先下降后又上升。而奎85-174在下降了10%時最高。在測定的時期中,該基因在石遠(yuǎn)321中表達(dá)量始終高于奎85-174。
GhGMPPA同樣是隨著脅迫程度的增加,呈現(xiàn)一個上調(diào)表達(dá)的趨勢,在石遠(yuǎn)321中該基因在土壤含水量下降了40% 時表達(dá)量最高。而奎85-174在下降了10%時 最高,兩個階段中抗旱性強(qiáng)的材料中的表達(dá)量顯著高于抗性弱的材料。GhPPT1在親本中表達(dá)情況是隨著時間的增加都有一個上調(diào)表達(dá)的趨勢,在石遠(yuǎn)321中該基因在土壤含水量下降了10% 時表達(dá)量最高,隨著脅迫程度的增加,表達(dá)量呈下降趨勢。而奎85-174在下降了10%時最高。而GhASP1在受到干旱脅迫后,隨這脅迫程度增加時,該基因呈現(xiàn)一個上調(diào)表達(dá)的趨勢,石遠(yuǎn)321中該基因在土壤含水量下降了50% 時表達(dá)量最高,隨著脅迫程度的增加,表達(dá)量呈上升趨勢。而在奎85-174中則在下降了10%時最高。GhRBCS與GhGMPPA相同,同樣是在遠(yuǎn)321中該基因在土壤含水量下降了40% 時表達(dá)量最高。而奎85-174在下降了10%時最高。并且在所有的測定時期中,5個基因在石遠(yuǎn)321中表達(dá)量始終高于奎85-174。
在受到水分脅迫時,抗旱性強(qiáng)的材料中這些基因的表達(dá)量均有上調(diào)表達(dá)的趨勢,并且抗旱性強(qiáng)的材料在輕度和中度水分脅迫的條件下的5個基因表達(dá)量顯著高于抗旱性弱的材料。圖6
注(Note):a:GhalkB,b:GhGMPPA,c:GhPPT1,d:GhASP1,e:GhRBCS
在干旱脅迫第3 d時,GhalkB在抗旱性強(qiáng)的材料中表達(dá)量顯著高于抗旱性弱的材料,并且在此時GhalkB的基因表達(dá)量達(dá)到最大值。由于呈快速上調(diào)表達(dá)的趨勢,該基因?qū)Ω珊得{迫前期較為敏感。而在抗旱性弱的材料中表達(dá)量變化不顯著。GhalkB在干旱脅迫前期起著重要作用,抗旱性強(qiáng)的材料在輕度脅迫時表現(xiàn)高表達(dá),抗旱性弱的材料表達(dá)量變化不顯著,該基因可能與棉花抗旱性有關(guān)。脅迫第3 d時,GhGMPPA在抗旱性強(qiáng)的材料中表達(dá)量顯著高于抗旱性弱的材料,也在此時達(dá)到最大值。受到脅迫后GhGMPPA基因呈現(xiàn)一個快速上調(diào)的表達(dá)趨勢,但是隨著脅迫程度的增加,呈現(xiàn)一個逐步下降的過程,而在抗旱性弱的材料中表達(dá)量的變化不顯著??购敌詮?qiáng)的材料中新陸早45號、貝爾斯諾、石遠(yuǎn)321三個材料在脅迫3 d時表達(dá)量達(dá)到最高,并且顯著高于抗旱性弱的材料,而ND359-5和KK153是在脅迫第6 d時表達(dá)量達(dá)到最大,并且顯著高于抗旱性弱的材料,該基因在不同材料中對干旱脅迫程度的敏感度不同;而在抗旱弱的材料中變化不顯著。同樣是在干旱脅迫第3 d(達(dá)到輕度脅迫條件)時,GhASP1基因表達(dá)量最高,并且顯著高于抗旱性弱的材料,也呈現(xiàn)快速上調(diào)表達(dá)的趨勢,在抗旱弱的材料中表達(dá)量變化不顯著。在干旱脅迫第6 d時(中度水分脅迫),在抗旱性材料中表達(dá)量達(dá)到最大,并且顯著高于抗旱性弱的材料,其中新陸早45號、貝爾斯諾在脅迫第3 d時的表達(dá)量也高于其它材料,呈現(xiàn)一個逐步上升的趨勢,該基因在不同材料中對干旱脅迫程度的敏感度不同;而在抗旱弱的材料中變化不顯著。圖7
注(Note):a:GhalkB,b:GhGMPPA,c:GhPPT1,d:GhASP1,e:GhRBCS
3.1Schneeberger等[5]利用混池測序在擬南芥中成功定位到導(dǎo)致葉片淺綠色突變體的基因。同時在擬南芥中還通過混池測序定位到配子突變導(dǎo)致不育的基因[6],以及突變導(dǎo)致延遲的基因[7]。Han等[29]在高粱莖含水量性狀相關(guān)研究匯總通過BSA測序分析定位關(guān)聯(lián)到6號染色體上339kb的區(qū)間內(nèi)[30]。賈秀萍等[31]通過高耐鹽自交系Y07-136R(父)和耐鹽敏感材聊Y05-222A(母)構(gòu)建的F2群體,通過BSA測序共發(fā)掘到6個耐鹽關(guān)鍵候選基因。劉夢雨等[32]用多油泡的羅紋金柑和少油泡的滑皮金柑雜交構(gòu)建F1群體,通過重測序-BSA分析共獲得了11個相關(guān)的候選基因。張堯峰等[33]在油菜中以有限花序突變株系FM8與野生型自交系 FM7雜交,F2中挑選子代構(gòu)建混池,獲得3個候選基因。研究中,以石遠(yuǎn)321為父本、奎85-174為母本構(gòu)建了陸地棉重組自交系。分別以重組自交系群體子代中,20份抗旱性強(qiáng)的子代材料和20份抗旱性弱的子代群體材料構(gòu)建混池,通過BSA測序分析,在SNP-index和ED兩種算法進(jìn)行關(guān)聯(lián)分析后取交集,最終關(guān)聯(lián)到3個相關(guān)的候選區(qū)域,候選區(qū)域的總長度為6.12Mb,其中共包含86個基因。結(jié)合GO注釋和KEGG通路分析,該基因在分子功能、細(xì)胞組分、或生物過程中有參與到干旱調(diào)控的過程中,參與了棉花的抗旱響應(yīng)過程。
3.2研究通過抗旱性評價篩選出抗旱性極端材料,結(jié)合BSA測序分析,篩選抗旱相關(guān)候選基因。通過測序數(shù)據(jù)分析及生物信息學(xué)分析,結(jié)合KEGG通路分析及基因注釋相關(guān)信息,共篩選出15個抗旱相關(guān)的候選基因,親本以及極端資源材料中熒光定量的測定中顯示,在不同程度的脅迫中,5個基因在親本和極端材料中表現(xiàn)出不同差異,其中alkB、GMPPA、PPT1、ASP1、RBCS基因的表達(dá)量在抗性強(qiáng)的材料中顯著高于抗旱性弱的材料,并且在前期抗旱性評價篩選出的5個關(guān)鍵指標(biāo)性狀的相關(guān)通路中,確定alkB、GMPPA、PPT1、ASP1、RBCS5個候選基因。通過數(shù)據(jù)庫比對分析,alkB是一種雙加氧酶通過影響水分運輸,導(dǎo)致葉片失水影響葉片的生長發(fā)育;GMPPA是GMPPB的同源物,可催化GDP-甘露糖的形成,后者是糖蛋白和糖脂的糖基部分的重要前體;PPT1作為磷酸烯醇丙酮酸轉(zhuǎn)運酶,將磷酸烯醇丙酮酸作為草酰乙酸途徑的一種底物導(dǎo)入葉綠體基質(zhì)中,同時參與葉肉細(xì)胞發(fā)育;ASP1作為氨基酸氨基轉(zhuǎn)移酶,參與氨基酸和三羧酸循環(huán),同時參與氮代謝及碳代謝方面;RBCS主要是在二氧化碳固定中將戊糖等進(jìn)行氧化裂解,參與光合作用。
得到3個關(guān)聯(lián)區(qū)域,總長度為6.12Mb,共包含86個基因。結(jié)合qRT-PCR分析,共篩選出GhalkB,GhGMPPA,GhPPT1,GhASP1,GhRBCS5個基因,以石遠(yuǎn)321正常對照cDNA為材料克隆了5個基因編碼序列,除GhPPT1存在跨膜區(qū)域,屬于膜蛋白外,其余基因均不存在跨膜區(qū)域,因此不屬于膜蛋白;其中GhRBCS亞細(xì)胞定位預(yù)測結(jié)果顯示定位于葉綠體中,預(yù)測定位位于線粒體中,而GhGMPPA和GhalkB則定位預(yù)測在細(xì)胞質(zhì)中,而GhPPT1則預(yù)測定位于葉綠體膜上;GhalkB基因含有Fe(2+)2-氧戊二酸雙加氧酶結(jié)構(gòu)域,與亞洲棉親緣關(guān)系最近;GhGMPPA基因在N末端域是GDP-甘露糖焦磷酸化酶同工型結(jié)構(gòu)域,其與雷蒙德氏棉親緣關(guān)系最近;GhPPT1基因具有磷酸三糖轉(zhuǎn)運蛋白家族蛋白結(jié)構(gòu)域,也與亞洲棉親緣關(guān)系最近;GhASP1基因具有PLP依賴性酶的天冬氨酸轉(zhuǎn)氨酶(AAT)結(jié)構(gòu)域,也與亞洲棉親緣關(guān)系最近;GhRBCS基因含有APF序列的保守結(jié)構(gòu)域,其中親緣關(guān)系最近的是海島棉和陸地棉的。在根、莖、葉中均有表達(dá),但是該5個基因在葉片中顯著表達(dá)。