基于BSA-seq 的黃褐棉纖維長(zhǎng)度候選基因發(fā)掘

陳 奇，王議萍，陳 艷，周婉婷，湯 備，朱心慧，于梓璇，李馨蕊，汪保華

（南通大學(xué) 生命科學(xué)學(xué)院，江蘇 南通 226019）

棉花是世界上重要的經(jīng)濟(jì)作物，也是最重要的天然纖維作物和紡織工業(yè)原料。纖維長(zhǎng)度是評(píng)價(jià)纖維品質(zhì)的重要指標(biāo)，在一定范圍內(nèi)，紗線強(qiáng)度隨著纖維長(zhǎng)度增加而增大；此外，相較于短纖維，長(zhǎng)纖維的加工效率更高，并能生產(chǎn)出質(zhì)量更高的細(xì)紗線[1]。唐淑榮等[2-3]研究指出，我國(guó)的栽培棉品種還存在一些突出的問(wèn)題，包括：斷裂比強(qiáng)度普遍偏低，纖維長(zhǎng)度、強(qiáng)度和細(xì)度搭配不合理，不能滿足棉紡織工業(yè)紡紗工藝的用棉需求，我國(guó)栽培棉的纖維品質(zhì)需進(jìn)一步提升。

棉花纖維長(zhǎng)度是由多基因控制的數(shù)量性狀，因此，定位其在染色體上的具體區(qū)域?qū)罄m(xù)候選基因的進(jìn)一步發(fā)掘具有重要意義。經(jīng)過(guò)多年研究，國(guó)內(nèi)外學(xué)者鑒定了很多棉花纖維長(zhǎng)度數(shù)量性狀基因座（quantitative trait locus，QTL），例如：在棉花1 號(hào)染色體上發(fā)現(xiàn)的一個(gè)穩(wěn)定的纖維長(zhǎng)度QTL，在不同研究團(tuán)隊(duì)的不同遺傳背景下都能檢測(cè)到。Chee 等[4]利用海島棉與陸地棉種間回交高世代群體定位到Chr.1染色體上的纖維長(zhǎng)度QTL（qFL-chr1），海島棉等位基因在3 個(gè)不同的回交高世代家系衍生的群體中均可提高纖維長(zhǎng)度，解釋12%～24%的表型方差。Shen 等[5]利用BC3F2植株構(gòu)建了3 個(gè)獨(dú)立的近等導(dǎo)入系群體，進(jìn)一步驗(yàn)證了該QTL；其中一個(gè)近等基因系R01-40-08 具有94.3%的輪回親本基因組，能顯著提高纖維長(zhǎng)度。Xu 等[6]利用R01-40-08 構(gòu)建了目標(biāo)區(qū)間分離的大群體，精細(xì)定位該QTL 至0.9 cM區(qū)間，并篩選到兩個(gè)候選基因GOBAR07705、GOBAR25992。除Chr.1 染色體上的穩(wěn)定主效QTL 外，研究人員也多次在Chr.7 染色體[7-9]和Chr.25 染色體[10-12]上檢測(cè)到纖維長(zhǎng)度QTL。此外，Ning 等[13]檢測(cè)到Chr.21 染色體上一個(gè)穩(wěn)定的纖維長(zhǎng)度QTL（qFLD11-1），在7 個(gè)環(huán)境下均能檢測(cè)到。Wang 等[14]利用一套陸地棉重組自交系群體鑒定到13 個(gè)纖維長(zhǎng)度QTL；其中，在Chr.18 染色體上檢測(cè)到纖維長(zhǎng)度QTL（qFL18.1）位于分子標(biāo)記PGML1145～JESPR178區(qū)間。雖然目前已鑒定了很多纖維長(zhǎng)度QTL，但克隆候選基因進(jìn)而得到應(yīng)用的并不多，因此，有必要進(jìn)一步發(fā)掘纖維長(zhǎng)度候選基因并驗(yàn)證其功能，最終用于棉花分子育種實(shí)踐。

棉屬分為A、B、C、D、E、F、G、K 8 個(gè)染色體組，共52 個(gè)種，其中包括45 個(gè)二倍體和7 個(gè)異源四倍體棉種[15]。在野生棉資源中，四倍體野生棉與陸地棉倍性相同，便于直接雜交獲得穩(wěn)定遺傳的后代；前人研究利用陸地棉與毛棉、陸地棉與達(dá)爾文棉種間群體，都鑒定到一些穩(wěn)定的纖維長(zhǎng)度QTL[16-17]。在棉花7 個(gè)異源四倍體里，黃褐棉與陸地棉遺傳距離最遠(yuǎn)[18]，因此，最有可能存在新的基因位點(diǎn)，挖掘并利用這些新的基因?qū)ν貙掙懙孛捱z傳背景和改良纖維品質(zhì)具有重要的科學(xué)意義和應(yīng)用價(jià)值。例如：Shen等[19]將陸地棉與黃褐棉雜交并回交，構(gòu)建了包含71個(gè)家系的染色體片段代換系群體，定位到29 個(gè)纖維品質(zhì)QTL，并發(fā)現(xiàn)部分QTL 的有利基因來(lái)自于黃褐棉。

本團(tuán)隊(duì)前期研究利用黃褐棉與陸地棉種間F2群體構(gòu)建了遺傳圖譜，在此基礎(chǔ)上構(gòu)建了高世代回交群體BC3F2、BC3F2:3和BC3F2:4，共鑒定出131 個(gè)控制纖維品質(zhì)性狀的QTL[20-23]?；诨亟桓呤来后w選擇了一套導(dǎo)入系，其中一些導(dǎo)入系纖維品質(zhì)表現(xiàn)突出；對(duì)這套導(dǎo)入系開(kāi)展了多年重復(fù)的田間試驗(yàn)，對(duì)纖維品質(zhì)性狀開(kāi)展了QTL 定位，共定位到15 個(gè)纖維長(zhǎng)度QTL，其中包括Chr.18 染色體上一個(gè)纖維長(zhǎng)度主效QTL（qUHM-18-1），其增效基因來(lái)自黃褐棉，解釋23.36%的表型方差，加性效應(yīng)為1.42 mm[24]，并由此通過(guò)分子標(biāo)記輔助選育出一個(gè)帶有qUHM-18-1且纖維長(zhǎng)度突出的黃褐棉導(dǎo)入系IL9（纖維長(zhǎng)度3 年平均值為33.3 mm）。本研究擬利用IL9 與其陸地棉輪回親本PD94042 構(gòu)建的F2大群體，通過(guò)極端性狀混池測(cè)序分析（bulked-segregant analysis sequencing，BSA-seq）結(jié)合QTL 定位結(jié)果，進(jìn)一步挖掘黃褐棉纖維長(zhǎng)度候選功能基因，為棉花分子育種提供基因資源。

1 實(shí)驗(yàn)材料與方法

1.1 植物材料

本實(shí)驗(yàn)室前期分子輔助選擇了一個(gè)優(yōu)質(zhì)黃褐棉導(dǎo)入系IL9，在IL9 的基因組構(gòu)成中，黃褐棉導(dǎo)入片段約占15.3%[24]。利用IL9 和陸地棉輪回親本PD94042 雜交并自交，構(gòu)建得到F2代次級(jí)群體，共1 244 個(gè)單株，開(kāi)展纖維品質(zhì)表型鑒定。

1.2 BSA-seq 分析

1.2.1 材料準(zhǔn)備

盛花期按單株采集BC3F2群體的棉花葉片放置于-80 ℃超低溫冰箱保存，在獲得纖維品質(zhì)表型數(shù)據(jù)后，根據(jù)測(cè)序要求挑選纖維長(zhǎng)度最長(zhǎng)和最短的棉花材料的葉片各30 份，再選取相同數(shù)量的親本材料PD94042 和IL9 的葉片，提取基因組DNA 組成4個(gè)混池，由北京組學(xué)生物科技有限公司進(jìn)行基因組重測(cè)序，測(cè)序深度為30×，參考基因組為棉花全基因組序列（http://mascotton.njau.edu.cn/Data.htm）。

1.2.2 文庫(kù)構(gòu)建及測(cè)序

首先，檢測(cè)樣品基因組的DNA，合格的DNA 用超聲波打斷的方法將其片段化；然后，進(jìn)行片段的純化和末端修復(fù)，進(jìn)而開(kāi)展3′端加A 和測(cè)序接頭的連接；之后，利用瓊脂糖凝膠電泳選擇片段大小，進(jìn)行PCR 擴(kuò)增形成測(cè)序文庫(kù)。構(gòu)建好的文庫(kù)先進(jìn)行質(zhì)檢，再利用Illumina HiSeq 平臺(tái)對(duì)質(zhì)檢合格的文庫(kù)進(jìn)行測(cè)序。

1.2.3 測(cè)序結(jié)果質(zhì)控

高通量測(cè)序得到的原始測(cè)序序列（Raw data 或Raw reads）包含序列信息及其對(duì)應(yīng)的測(cè)序質(zhì)量信息，必須經(jīng)過(guò)質(zhì)控進(jìn)而得到Clean data 后才能用于后續(xù)分析。質(zhì)控方法包括：低質(zhì)量數(shù)據(jù)過(guò)濾、堿基測(cè)序質(zhì)量分布和堿基類型分布。

1.2.4 候選基因的發(fā)掘與鑒定

將質(zhì)控后得到的Clean data 重新比對(duì)到參考基因組上，進(jìn)而開(kāi)展后續(xù)的變異分析，根據(jù)比對(duì)結(jié)果進(jìn)行插入缺失（insertion-deletion，InDel）和單核苷酸多態(tài)性（single nucleotide polymorphism，SNP）的檢測(cè)及注釋。對(duì)于InDel 和SNP 的檢測(cè)主要使用GATK軟件工具包來(lái)實(shí)現(xiàn)，主要檢測(cè)過(guò)程可參考GATK 官方網(wǎng)站的BestPractice（https://www.broadinstitute.org/gatk/guide/best-practices.php），隨后使用一種高效的軟件工具ANNOVAR 對(duì)SNP 進(jìn)行注釋。對(duì)于差異基因的關(guān)聯(lián)分析，本實(shí)驗(yàn)通過(guò)兩種方法來(lái)確定候選區(qū)域，分別為歐式距離（euclidean distance，ED）算法關(guān)聯(lián)分析和基于BSA-seq 的SNP-index 算法關(guān)聯(lián)分析。SNP 質(zhì)量分布和SNP 突變頻譜被用來(lái)檢測(cè)所得SNP 數(shù)據(jù)的可靠性，最終根據(jù)SNP-index 算法計(jì)算與性狀關(guān)聯(lián)的候選區(qū)域。SNP-index 算法如下：基于基因分型的結(jié)果，篩選親本中純合的且親本間具有多態(tài)性的SNP 位點(diǎn)。將其中的一個(gè)親本作為參考，計(jì)算兩個(gè)子代在親本間標(biāo)記位點(diǎn)的SNP-index（SNP頻率），即統(tǒng)計(jì)混池和親本（或參考基因組）在某一個(gè)堿基位點(diǎn)上相同或不同的reads 數(shù)，進(jìn)而計(jì)算出不同reads 數(shù)占總reads 數(shù)的比例，即為該堿基位點(diǎn)的SNP-index。完全與其相同的SNP-index 為0；完全與其不同的SNP-index 為1。計(jì)算Δ（SNP-index），即兩個(gè)子代的SNP-index 之差：Δ（SNP-index）=SNP-index（極端性狀B）-SNP-index（極端性狀A(yù)）。該差值越接近于1，則標(biāo)記SNP 與目標(biāo)性狀之間的關(guān)聯(lián)度越大，把超過(guò)閾值線的區(qū)域定為候選區(qū)域。對(duì)分析所得的候選基因區(qū)域內(nèi)的基因進(jìn)行基因同源性分析及相關(guān)功能深度注釋，以此來(lái)篩選所要得到的候選基因。

2 結(jié)果與分析

2.1 表型結(jié)果統(tǒng)計(jì)分析

利用本實(shí)驗(yàn)室前期分子輔助選擇的優(yōu)質(zhì)黃褐棉導(dǎo)入系IL9 和陸地棉輪回親本PD94042 雜交構(gòu)建得到F2代次級(jí)群體，共1 244 個(gè)單株，在相同條件下，對(duì)親本材料PD94042、IL9 和F2代群體的纖維長(zhǎng)度表型數(shù)據(jù)進(jìn)行測(cè)量和記錄。該部分表型數(shù)據(jù)被用來(lái)進(jìn)行BSA 分析所需數(shù)據(jù)的選擇，F(xiàn)2代次級(jí)群體平均表型結(jié)果如表1 所示。從表中的數(shù)據(jù)可以看出，相較于親本PD94042，導(dǎo)入系IL9 明顯具有更優(yōu)良的纖維表型表現(xiàn)，后代群體的纖維長(zhǎng)度得到了改良。

表1 F2 代次級(jí)群體纖維長(zhǎng)度表型數(shù)據(jù)Tab.1 Phenotypic data of fiber length in the F2 population

2.2 與參考基因組的比對(duì)統(tǒng)計(jì)

將重測(cè)序獲得的reads 重新比對(duì)到參考基因組上，進(jìn)而開(kāi)展后續(xù)的變異分析。本研究的測(cè)序樣品序列reads 總數(shù)在277～327 Mbp 之間，GC 堿基數(shù)都在35%左右，GC 分布正常，而且各個(gè)樣品Q30 的值都在90%左右，說(shuō)明過(guò)濾后的結(jié)果較為可靠。

經(jīng)過(guò)過(guò)濾和修飾讀取后，4 個(gè)樣本reads 比對(duì)到參考基因組上的覆蓋率都超過(guò)98%，雙端測(cè)序序列均比對(duì)到參考基因組上，且距離符合測(cè)序片段的長(zhǎng)度分布均超過(guò)94%。對(duì)于兩個(gè)極端性狀DNA 混池，測(cè)序深度達(dá)到33×，且參考基因組中至少有10 個(gè)堿基覆蓋的位點(diǎn)占基因組的百分比均超過(guò)81%，每個(gè)樣品的測(cè)序數(shù)據(jù)信息充分，測(cè)序質(zhì)量合格。結(jié)果表明該數(shù)據(jù)符合突變分析檢測(cè)的標(biāo)準(zhǔn)，可成功用于Illumina 測(cè)序分析。

2.3 混池SNP 檢測(cè)統(tǒng)計(jì)

樣品表型之間的差異往往是其基因組序列上的變異所導(dǎo)致的，在BSA-seq 分析中，利用GATK軟件工具包，把重心放在檢測(cè)樣品間序列的SNP上。SNP 指由基因組核苷酸水平上的變異引起的DNA 序列多態(tài)性，包括由單堿基的轉(zhuǎn)換、顛換和單堿基的插入與缺失等引起的同義突變和錯(cuò)義突變。通過(guò)與參考基因組的對(duì)比，得到每一個(gè)樣本和參考基因組的SNP 數(shù)目，以及樣本之間的結(jié)果比較。混池間的SNP 結(jié)果統(tǒng)計(jì)如表2 所示。

表2 混池SNP 結(jié)果統(tǒng)計(jì)表Tab.2 SNP information of different samples

從表2 中可以看出，F(xiàn)Lmax中的SNP 總數(shù)高達(dá)2 829 100，而FLmin中也包含1 885 065 個(gè)SNP。為進(jìn)一步統(tǒng)計(jì)和可視化4 個(gè)樣本間的SNP 重復(fù)性，我們做了SNP 的統(tǒng)計(jì)結(jié)果Venn 圖（圖1）。從Venn 圖整體分布情況來(lái)看，F(xiàn)Lmax混池基因組更接近于IL9，F(xiàn)Lmin混池基因序列則更接近于PD94042。這也從側(cè)面說(shuō)明了子代中控制更長(zhǎng)纖維長(zhǎng)度的基因大部分來(lái)源于黃褐棉導(dǎo)入系IL9，該結(jié)果為本研究課題的展開(kāi)和后續(xù)工作提供了理論支持。

圖1 纖維長(zhǎng)度混池樣品SNP 統(tǒng)計(jì)Venn 圖Fig.1 Venn diagram of fiber length SNPs

2.4 候選基因區(qū)域篩選

本研究通過(guò)基于BSA-seq 的SNP-index 關(guān)聯(lián)分析來(lái)計(jì)算兩個(gè)混池間等位基因的基因型頻率，進(jìn)而獲得與目標(biāo)性狀相關(guān)聯(lián)的基因組區(qū)域。根據(jù)SNP-index 方法關(guān)聯(lián)閾值判定，首先計(jì)算兩個(gè)子代在親本間標(biāo)記位點(diǎn)的SNP 的頻率（即SNP-index），計(jì)算Δ（SNP-index），即兩個(gè)子代的SNP-index 作差，超過(guò)閾值線的為候選基因組區(qū)域。當(dāng)擬合后的Δ（SNP-index）置信度為0.99 時(shí)，共得到3 個(gè)與棉花纖維長(zhǎng)度性狀相關(guān)聯(lián)的染色體區(qū)域，分布于Chr.7、Chr.10、Chr.18 3 條染色體上（圖2）。

圖2 纖維長(zhǎng)度混池SNP-index 關(guān)聯(lián)分析Fig.2 SNP-index association analysis of fiber length

通過(guò)對(duì)關(guān)聯(lián)分析所得到的3 個(gè)候選區(qū)域的進(jìn)一步統(tǒng)計(jì)，包括這些區(qū)域所處的染色體位置（起始和結(jié)束位點(diǎn)）、序列長(zhǎng)度和候選基因數(shù)目，得到候選區(qū)域的相關(guān)信息統(tǒng)計(jì)（表3）。從表中數(shù)據(jù)可以看出，纖維長(zhǎng)度性狀候選基因主要集中在Chr.7、Chr.10、Chr.18 3 條染色體的3 個(gè)基因區(qū)域，涉及的基因數(shù)分別為7、21、6 個(gè)。

表3 纖維長(zhǎng)度基因候選區(qū)域相關(guān)信息統(tǒng)計(jì)表Tab.3 Statistics of candidate gene regions related to fiber length

2.5 候選基因篩選

在前期研究中，我們鑒定了一個(gè)位于18 號(hào)染色體上的纖維長(zhǎng)度主效QTL（qUHM-18-1），候選區(qū)域覆蓋染色體的31 696 068～51 883 215 bp 基因組區(qū)間，而本次BSA-seq 研究中鑒定到的Chr.18（D13）候選基因區(qū)域與該QTL 的染色體區(qū)域有部分重疊，且重疊區(qū)域包含了引起氨基酸變化的單核苷酸變化，因此我們選取Chr.18 染色體作為后續(xù)驗(yàn)證的主要候選基因區(qū)域。我們對(duì)BSA 候選區(qū)域內(nèi)的基因開(kāi)展了進(jìn)一步的整理，隨后通過(guò)GO 功能注釋信息和KEGG 通路注釋信息，進(jìn)而對(duì)區(qū)域內(nèi)的候選基因進(jìn)行與纖維長(zhǎng)度性狀相關(guān)的篩選。綜合QTL 及BSA-seq 結(jié)果，該區(qū)域存在6 個(gè)相關(guān)基因，分別為Gh_D13G1292、Gh_D13G1293、Gh_D13G1294、Gh_D 13G1295、Gh_D13G1296、Gh_D13G1297。通過(guò)GO功能注釋信息和KEGG 通路分析，最終篩選到3 個(gè)候選基因：Gh_D13G1294、Gh_D13G1295 和Gh_D 13G1296（基因序列信息見(jiàn)表4）。

表4 目的基因序列信息Tab.4 Sequences of the target genes

3 結(jié)論

前期研究中，我們基于一套黃褐棉導(dǎo)入系群體開(kāi)展了纖維品質(zhì)QTL 定位，鑒定了15 個(gè)纖維長(zhǎng)度主效QTL[24]。選擇纖維品質(zhì)突出的導(dǎo)入系IL9 與其陸地棉輪回親本PD94042 雜交并自交，構(gòu)建了F2大群體，結(jié)合表型鑒定結(jié)果，開(kāi)展了纖維長(zhǎng)度的BSA-seq 分析?；赟NP-index 和ED 的關(guān)聯(lián)算法，控制纖維長(zhǎng)度的數(shù)量性狀位點(diǎn)被分別定位在Chr.7、Chr.10、Chr.18 共3 條染色體上3 個(gè)相關(guān)區(qū)域內(nèi)，涉及的基因數(shù)分別為7、21、6 個(gè)。前期基于黃褐棉導(dǎo)入系的QTL 定位研究鑒定到Chr.18 上一個(gè)纖維長(zhǎng)度QTL（qUHM-18-1），其增效基因來(lái)自黃褐棉，解釋23.36%的表型方差[24]。該QTL 是在前人研究發(fā)現(xiàn)的來(lái)自達(dá)爾文棉、毛棉等四倍體野生棉種之外，增效基因來(lái)自黃褐棉的新的纖維長(zhǎng)度主效QTL，為陸地棉纖維品質(zhì)改良提供了新的基因資源。通過(guò)比較該QTL 和BSA-seq 候選區(qū)域的位置信息，我們將候選基因區(qū)域縮小到更精細(xì)的區(qū)域。對(duì)分析所得區(qū)域內(nèi)的候選基因進(jìn)行基因同源性分析及相關(guān)功能注釋篩選，最終篩選出最可能和纖維長(zhǎng)度的生長(zhǎng)發(fā)育建成相關(guān)的3 個(gè)候選基因，為Gh_D13G1294、Gh_D 13G1295 和Gh_D13G1296。

Gh_D13G1294 編碼MYB-like102（MYB102），在擬南芥中，MYB102 轉(zhuǎn)錄因子可以通過(guò)調(diào)節(jié)乙烯生物合成[25]和細(xì)胞分裂素水平[26]來(lái)促進(jìn)植物體的生長(zhǎng)發(fā)育，但在棉花中的分子作用機(jī)制尚不明確。Gh_D13G1295 編碼myb domain protein 4r1 轉(zhuǎn)錄因子，屬于R2R3-MYB 轉(zhuǎn)錄因子基因家族，是植物中最大的轉(zhuǎn)錄因子基因家族之一，該基因家族在植物的生長(zhǎng)發(fā)育建成和代謝調(diào)節(jié)方面發(fā)揮廣泛作用。Gh_D13G1296 編 碼myb domain protein 16 轉(zhuǎn)錄因子（MYB16），該轉(zhuǎn)錄因子在植物體的表皮形成[27]和細(xì)胞形態(tài)變化發(fā)生[28]中發(fā)揮著重要作用。在分子功能方面，這三個(gè)基因都是蛋白質(zhì)編碼基因，且都與染色質(zhì)結(jié)合（GO:0003682）相關(guān)，通過(guò)編碼相關(guān)蛋白與染色質(zhì)、蛋白質(zhì)、DNA 的纖維網(wǎng)絡(luò)或者在染色體間期組成真核細(xì)胞核染色體的RNA 進(jìn)行選擇性和非共價(jià)相互作用完成相關(guān)功能。KEGG 通路分析發(fā)現(xiàn)這幾個(gè)基因都主要集中作用在DNA 的結(jié)合過(guò)程中，以及作用于其合成蛋白質(zhì)的轉(zhuǎn)錄因子方面。這些基因可能在棉花纖維長(zhǎng)度的發(fā)育過(guò)程中起重要作用，在后續(xù)研究中，將進(jìn)一步驗(yàn)證其功能，最終用于棉花分子育種實(shí)踐。