基于BSA分析技術對陸地棉抗旱基因的鑒定

孫豐磊,任姣姣,雷 斌,高文偉,曲延英

(1.新疆農業大學農學院,烏魯木齊 830052;2 新疆農業科學院核技術生物研究所,烏魯木齊 830052)

0 引 言

【研究意義】不同棉花品種的不同生長發育階段對干旱脅迫響應不同。棉花花鈴期若受到干旱脅迫會影響棉花產量。以現有棉花材料為基礎,篩選抗旱性強材料,挖掘、研究棉花抗旱的相關基因,對提高棉花的抗旱性及加快培育抗旱棉花新品種的進程有重要意義。【前人研究進展】作物抗旱性是一個復雜的數量性狀,傳統的定位方法耗時長、成本高[1]。混池測序是一種基于二代測序以集群分離方法為主的技術,應用于基因定位和克隆等[2]。BSA(集群分離方法)應用在萵苣的抗霜霉病研究中,篩選到3個RAPD分子標記[3],3個分子標記與抗性基因Dm5/8緊密連鎖[4]。通過選擇具有極端性狀或具有明顯表型差異的親本來構建群體,從子代群體中根據目標性狀的特點選擇具有極端差異子代個體,構成兩個極端性狀的子代DNA混池,通過測序分析獲得與性狀相關的分子標記,定位相關基因。由于作為兩個親本是在目標性狀中選擇明顯差異的個體,因此,在子代的兩個極端混池中,只有在目標性狀基因的區域中有差異,其它區域的差異較小或基本沒有差異,增加了定位的準確性和可靠性。混池測序(BSA-seq)早期應用在基因組較小的物種中,并且目標基因在該物種的突變體中。Schneeberger等[5]利用混池測序成功定位到導致擬南芥葉片淺綠色突變體的基因。在擬南芥中還定位到配子突變導致不育的基因[6]以及突變導致開花延遲的基因[7]。而此方法在水稻中也得到應用,Abe等[8]在EMS誘變的突變體水稻中將白化和半矮稈定位在某一區間內,進而對目標進行了定位和克隆。Takagi等[9]在EMS誘變的水稻突變體中,定位到一個耐鹽堿基因,通過MutMap 的方法克隆該基因并培育出耐鹽堿的品種。陳忠明等[10]結合BSA分析法和RFLP的方法,在可育與半不育的極端混池中發現了2個標記(RG213、G200)與親和基因連鎖。陸光遠等[11]利用油菜的回交群體篩選到兩個AFLP標記,2個標記與細胞核雄性不育基因緊密連鎖。李傳友等[12]結合BSA和AFLP分析法篩選到光敏核不育基因相關的多態性產物。張顯等[13]在西瓜隱形和雄性不育和可育的植株中,通過RAPD-BSA方法篩選出s1167的引物可以區分該兩種株系。Wang等[14]通過BSA分析的方法在水稻中篩選到9個分子標記,利用9個分子標記最終定位到tms5基因。Geng等[15]利用BSA分析法對油菜千粒重基因進行定位,篩選到4個與千粒重相關的候選基因。Hong等[16]通過該方法篩選到與甜味相關的候選基因1個,酸味相關的候選基因8個。Zheng等[17]通過混池測序,共篩選到4個抗稻瘟病候選基因。而Zhang等[18]通過混池關聯分析的方法研究黃瓜抗白粉病相關基因,共篩選到6個相關的候選基因。Chen等[19]和Feng等[20]通過BSA測序分析,再結合SSR標記,定位到與抗稻瘟病相關的基因。Qin等[21]在大麥白化性狀方面利用集群分離分析的方法,結合SSR分子標記,定位6個有關該性狀的候選基因。杜威世等[22]在海島棉品種和陸地棉品種的雜交F2代中,結合SSR分子標記和BSA法定位到3個與抗棉花黃萎病相關的基因標記。Hayes等[23]研究結果表明,在大豆中通過BSA法篩選到4個Rsv1相關的分子標記,該4個分子標記與抗花葉病相關。【本研究切入點】棉花抗旱性是一個復雜的數量性狀,傳統的定位方法耗時長、成本高。混池測序是一種基于二代測序以集群分離方法為主的技術,可以有效用于棉花復雜性狀的遺傳定位。【擬解決的關鍵問題】篩選出抗旱性強和抗旱性敏感的材料,構建抗旱性強和抗旱性弱的兩個極端樣本混池,采用BSA混池測序技術測序分析親本及子代混池,定位篩選抗旱基因位點候選區域,為挖掘棉花相關的抗旱基因奠定基礎。

1 材料與方法

1.1 材 料

選用由石遠321為父本和奎85-174為母本構建的150份材料(重組自交系子代群體2年2點田間抗旱表型鑒定),篩選得到抗旱性強的20份極端材料和20份水分敏感的極端材料,分別提取極端材料的葉片DNA,檢測濃度以及純度,等量混合,分別構成極端抗旱池與極端敏旱池,提取2個親本的葉片DNA。

樣本混池建庫和測序由北京百邁客生物科技有限公司完成。按照lllumina公司執行標準,首先檢測DNA樣品后,進行超聲破碎[24],修復、拼接、擴增,并構建完成測序所需要的文庫。經過質檢合格后的文庫用測序分析(參考陸地棉基因組(版本號為v2.1))。

1.2 方 法

1.2.1 SNP變異位點的檢測

SNP變異位點檢測是通過GATK[25-26]軟件工具包完成。根據Clean Reads在參考基因組的定位結果,通過去重復(Mark Duplicates)以及GATK進行局部重比對(Local Realignment)、堿基質量值校正(Base Recalibration)等處理,再使用GATK檢測單核苷酸多態性(Single Nucleotide Polymorphism,SNP),得到SNP位點集。

根據Clean Reads在參考基因組上的定位結果,小片段的插入與缺失有基因組之間的檢測確定。但同樣反映了其差異的Small InDel變異的數量較少,其中InDel時會引起移碼突變[27],如果移碼突變區在編碼區將導致基因功能上的改變。

1.2.2 關聯分析及候選區域的基因注釋

對SNP過濾后,進行關聯分析。在SNP關聯分析中主要采用歐式距離(Euclidean Distance,ED)和SNP-inde 2種分析方法,將2種分析方法共同關聯到的區間作為抗旱相關的主要區間。

其中歐氏距離分析方法(Euclidean Distance,ED)是通過分析測序數據在混池間確定差異顯著的標記,確定與目標性狀關聯的區域。在目標性狀相關位點上除了2個子代混池與親本間存在顯著差異外,其它位點不存在差異或趨于基本一致,通過分析混池間差異的SNP位點,計算ED值,在實驗中取ED的2次方進行擬合。

1.2.3 熒光定量驗證候選基因的差異表達

將10份抗旱性極端材料在實驗室內進行土培(土壤為黑土∶蛭石為2∶1混合)培養至三葉期時,進行土壤水分脅迫試驗,用于熒光定量測定基因表達量。

1.3 數據處理

統計混池間基因型頻率差異的分析方法應用在標記關聯分析中稱為SNP-index算法[28]。是通過Δ(SNP-index)統計分析混池間基因型頻率間的顯著差異。選擇性狀相關的區域是在標記的染色體上通過ΔSNP-index值進行擬合分析后關聯閾值以上的區域。

通過BLAST軟件對候選基因在多個數據中(GO[20]、KEGG[21]、NR[22])比對詳細注釋。

2 結果與分析

2.1 樣品測序數據對比

研究表明,測序得到原始reads為119 060 015、121 194 054、190 811 397、223 420 448,得到clean reads與原始reads基本一致,總數據量為197.65 Gbp,其中Q30(%)在93%以上,并且GC含量在35.43%～35.73%,共得到142.23 Gbp的Clean reads,平均每個樣品測序深度15.47X,最終的平均比對效率為99.045%,達到18.75X平均覆蓋深度和99.03%的基因組覆蓋度。表1,表2

表1 樣品測序數據質量統計

表2 與參考基因比對

2.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 抗旱性關聯

2.3.1 SNP檢測

研究表明,共得到2 561 910個SNP位點,再次過濾后最終得到高質量SNP位點共540 596個。

取所有位點的擬合值作為關聯閾值,得到ED值為0.17,共得到10個區域,總長度為39.30 Mb,共包含979個基因,其中包含非同義突變SNP位點的基因共112個。圖2,表3

注:橫坐標為染色體名稱,彩色的點代表每個SNP位點的ED值,黑色的線為擬合后的ED值,紅色的虛線代表顯著性關聯閾值,下同

表3 ED算法關聯候選區域

擬合后的ΔSNP-index值為0.31,共得到5個區域,總長度為6.12 Mb,共包含86個基因,在所有基因非同義突變的共有27個。圖3,表4

圖3 ΔSNP-index關聯值在染色體上分布

表4 SNP-index算法關聯候選區域

2.3.2 候選區域篩選及候選區域基因注釋

研究表明,共篩選到3個候選區域,總長度為6.12 Mb,共包含86個基因。表5

表5 ED和SNP-index兩種算法關聯候選區域

親本間存在非同義突變的SNP共43個,混池間存在非同義突變的SNP共20個,在候選區域內共注釋到83個基因,而注釋到的非同義突變基因在親本間共有27個。其中在GO數據庫中注釋到49個基因,而在KEGG數據庫中注釋到23個基因。

候選區域內每個組分中基因數量不同,其中大多數被劃分到生物過程(Biological process)中。前5的功能分類分別是:在細胞組分中為質體部分、葉綠體部分、細胞器膜、轉移酶復合體和葉綠體;在分子功能中為RNA導向的DNA聚合酶活性、核糖二磷酸羧化酶活性、蛋白質結構域特異性結合、雙加氧酶活性、碳水化合物磷酸酶活性;生物過程中主要為DNA整合、依賴RNA的DNA復制、糖醇代謝過程、細胞壁組織或生物合成、有機羥基化合物代謝過程。圖4,表6

圖4 候選區域內基因GO注釋聚類

表6 候選區域內基因GO注釋

得到的代謝通路中基因注釋分布最多的在代謝作用和遺傳信息處理中,代謝作用通路中注釋基因最多的是內質網中的蛋白質加工和RNA聚合酶;在遺傳信息處理的代謝通路上苯丙烷生物合成和氨基酸生物合成中注釋到的基因種類最多。圖5

2.3.3 候選基因的篩選

研究表明,參與細胞組成的相關基因與棉花的抗旱性相關,包括細胞核(GO:0005634),葉綠體(GO:0009507),過氧化物酶體(GO:0005777),RNA聚合酶Ⅲ復合物(GO:0005666)等,GH_A07G0705、GH_A07G0701、GH_D07G1930、GH_D07G1932、GH_D07G1934、GH_A07G0708、GH_A07G0695、GH_A07G0703等8個基因參與調控棉花的抗旱性。

鋅離子結合(GO:0008270)、RNA定向的DNA聚合酶活性(GO:0003964)、蛋白質結構域特異性結合(GO:0019904)、甲基轉移酶活性(GO:0008168)等參與到干旱脅迫相關反應中。GH_A07G0695、GH_A10G1602、GH_D07G1927、GH_A10G1623、GH_A10G1601等24個基因與干旱脅迫調控相關。

一些參與干旱脅迫調控的過程,如代謝過程(GO:0008152),響應脫落酸反應過程(GO:0009737),有機物的運輸過程(GO:0071702),還有氧化還原過程(GO:0055114)等均參與干旱脅迫的生物調控過程中。預測GH_A07G0696、GH_A07G0700、GH_A07G0695、GH_A07G0703、GH_A07G0695、GH_A07G0696、GH_A07G0705等24個基因與干旱脅迫相關。表7

篩選出與前期5個棉花抗旱性評價關鍵性狀指標相關的15個候選基因,分布在3個候選區域內。表8

表8 抗旱候選基因及注釋

2.4 qRT-PCR分析候選基因表模式

研究表明,在水分脅迫達到輕度(土壤含水下降了20%)、中度(土壤含水下降了40%)、重度(土壤含水下降了50%)時,在脅迫程度增加時,GhalkB隨著脅迫程度的增加都有一個上調表達的趨勢,在石遠321中該基因在土壤含水量下降了10% 時表達量最高,隨著脅迫程度的增加,該基因的表達量先下降后又上升。而奎85-174在下降了10%時最高。在測定的時期中,該基因在石遠321中表達量始終高于奎85-174。

GhGMPPA同樣是隨著脅迫程度的增加,呈現一個上調表達的趨勢,在石遠321中該基因在土壤含水量下降了40% 時表達量最高。而奎85-174在下降了10%時 最高,兩個階段中抗旱性強的材料中的表達量顯著高于抗性弱的材料。GhPPT1在親本中表達情況是隨著時間的增加都有一個上調表達的趨勢,在石遠321中該基因在土壤含水量下降了10% 時表達量最高,隨著脅迫程度的增加,表達量呈下降趨勢。而奎85-174在下降了10%時最高。而GhASP1在受到干旱脅迫后,隨這脅迫程度增加時,該基因呈現一個上調表達的趨勢,石遠321中該基因在土壤含水量下降了50% 時表達量最高,隨著脅迫程度的增加,表達量呈上升趨勢。而在奎85-174中則在下降了10%時最高。GhRBCS與GhGMPPA相同,同樣是在遠321中該基因在土壤含水量下降了40% 時表達量最高。而奎85-174在下降了10%時最高。并且在所有的測定時期中,5個基因在石遠321中表達量始終高于奎85-174。

在受到水分脅迫時,抗旱性強的材料中這些基因的表達量均有上調表達的趨勢,并且抗旱性強的材料在輕度和中度水分脅迫的條件下的5個基因表達量顯著高于抗旱性弱的材料。圖6

注(Note):a:GhalkB,b:GhGMPPA,c:GhPPT1,d:GhASP1,e:GhRBCS

在干旱脅迫第3 d時,GhalkB在抗旱性強的材料中表達量顯著高于抗旱性弱的材料,并且在此時GhalkB的基因表達量達到最大值。由于呈快速上調表達的趨勢,該基因對干旱脅迫前期較為敏感。而在抗旱性弱的材料中表達量變化不顯著。GhalkB在干旱脅迫前期起著重要作用,抗旱性強的材料在輕度脅迫時表現高表達,抗旱性弱的材料表達量變化不顯著,該基因可能與棉花抗旱性有關。脅迫第3 d時,GhGMPPA在抗旱性強的材料中表達量顯著高于抗旱性弱的材料,也在此時達到最大值。受到脅迫后GhGMPPA基因呈現一個快速上調的表達趨勢,但是隨著脅迫程度的增加,呈現一個逐步下降的過程,而在抗旱性弱的材料中表達量的變化不顯著。抗旱性強的材料中新陸早45號、貝爾斯諾、石遠321三個材料在脅迫3 d時表達量達到最高,并且顯著高于抗旱性弱的材料,而ND359-5和KK153是在脅迫第6 d時表達量達到最大,并且顯著高于抗旱性弱的材料,該基因在不同材料中對干旱脅迫程度的敏感度不同;而在抗旱弱的材料中變化不顯著。同樣是在干旱脅迫第3 d(達到輕度脅迫條件)時,GhASP1基因表達量最高,并且顯著高于抗旱性弱的材料,也呈現快速上調表達的趨勢,在抗旱弱的材料中表達量變化不顯著。在干旱脅迫第6 d時(中度水分脅迫),在抗旱性材料中表達量達到最大,并且顯著高于抗旱性弱的材料,其中新陸早45號、貝爾斯諾在脅迫第3 d時的表達量也高于其它材料,呈現一個逐步上升的趨勢,該基因在不同材料中對干旱脅迫程度的敏感度不同;而在抗旱弱的材料中變化不顯著。圖7

注(Note):a:GhalkB,b:GhGMPPA,c:GhPPT1,d:GhASP1,e:GhRBCS

3 討 論

3.1Schneeberger等[5]利用混池測序在擬南芥中成功定位到導致葉片淺綠色突變體的基因。同時在擬南芥中還通過混池測序定位到配子突變導致不育的基因[6],以及突變導致延遲的基因[7]。Han等[29]在高粱莖含水量性狀相關研究匯總通過BSA測序分析定位關聯到6號染色體上339kb的區間內[30]。賈秀萍等[31]通過高耐鹽自交系Y07-136R(父)和耐鹽敏感材聊Y05-222A(母)構建的F2群體,通過BSA測序共發掘到6個耐鹽關鍵候選基因。劉夢雨等[32]用多油泡的羅紋金柑和少油泡的滑皮金柑雜交構建F1群體,通過重測序-BSA分析共獲得了11個相關的候選基因。張堯峰等[33]在油菜中以有限花序突變株系FM8與野生型自交系 FM7雜交,F2中挑選子代構建混池,獲得3個候選基因。研究中,以石遠321為父本、奎85-174為母本構建了陸地棉重組自交系。分別以重組自交系群體子代中,20份抗旱性強的子代材料和20份抗旱性弱的子代群體材料構建混池,通過BSA測序分析,在SNP-index和ED兩種算法進行關聯分析后取交集,最終關聯到3個相關的候選區域,候選區域的總長度為6.12Mb,其中共包含86個基因。結合GO注釋和KEGG通路分析,該基因在分子功能、細胞組分、或生物過程中有參與到干旱調控的過程中,參與了棉花的抗旱響應過程。

3.2研究通過抗旱性評價篩選出抗旱性極端材料,結合BSA測序分析,篩選抗旱相關候選基因。通過測序數據分析及生物信息學分析,結合KEGG通路分析及基因注釋相關信息,共篩選出15個抗旱相關的候選基因,親本以及極端資源材料中熒光定量的測定中顯示,在不同程度的脅迫中,5個基因在親本和極端材料中表現出不同差異,其中alkB、GMPPA、PPT1、ASP1、RBCS基因的表達量在抗性強的材料中顯著高于抗旱性弱的材料,并且在前期抗旱性評價篩選出的5個關鍵指標性狀的相關通路中,確定alkB、GMPPA、PPT1、ASP1、RBCS5個候選基因。通過數據庫比對分析,alkB是一種雙加氧酶通過影響水分運輸,導致葉片失水影響葉片的生長發育;GMPPA是GMPPB的同源物,可催化GDP-甘露糖的形成,后者是糖蛋白和糖脂的糖基部分的重要前體;PPT1作為磷酸烯醇丙酮酸轉運酶,將磷酸烯醇丙酮酸作為草酰乙酸途徑的一種底物導入葉綠體基質中,同時參與葉肉細胞發育;ASP1作為氨基酸氨基轉移酶,參與氨基酸和三羧酸循環,同時參與氮代謝及碳代謝方面;RBCS主要是在二氧化碳固定中將戊糖等進行氧化裂解,參與光合作用。

4 結 論

得到3個關聯區域,總長度為6.12Mb,共包含86個基因。結合qRT-PCR分析,共篩選出GhalkB,GhGMPPA,GhPPT1,GhASP1,GhRBCS5個基因,以石遠321正常對照cDNA為材料克隆了5個基因編碼序列,除GhPPT1存在跨膜區域,屬于膜蛋白外,其余基因均不存在跨膜區域,因此不屬于膜蛋白;其中GhRBCS亞細胞定位預測結果顯示定位于葉綠體中,預測定位位于線粒體中,而GhGMPPA和GhalkB則定位預測在細胞質中,而GhPPT1則預測定位于葉綠體膜上;GhalkB基因含有Fe(2+)2-氧戊二酸雙加氧酶結構域,與亞洲棉親緣關系最近;GhGMPPA基因在N末端域是GDP-甘露糖焦磷酸化酶同工型結構域,其與雷蒙德氏棉親緣關系最近;GhPPT1基因具有磷酸三糖轉運蛋白家族蛋白結構域,也與亞洲棉親緣關系最近;GhASP1基因具有PLP依賴性酶的天冬氨酸轉氨酶(AAT)結構域,也與亞洲棉親緣關系最近;GhRBCS基因含有APF序列的保守結構域,其中親緣關系最近的是海島棉和陸地棉的。在根、莖、葉中均有表達,但是該5個基因在葉片中顯著表達。