基于BSA-seq的擬南芥輻射誘變突變體的基因定位

汪澤民 何 曦 張宇涵 周 明 洪麗蘭，*

（1浙江大學(xué)原子核農(nóng)業(yè)科學(xué)研究所/農(nóng)業(yè)農(nóng)村部和浙江省核農(nóng)學(xué)重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310058；2浙江大學(xué)生命科學(xué)學(xué)院植物生物學(xué)研究所，浙江 杭州 310058）

從人工誘導(dǎo)或自然發(fā)生的突變體中挖掘出突變信息既是分子遺傳研究的重要內(nèi)容，也是培育高產(chǎn)優(yōu)質(zhì)新品種的重要途徑［1-2］。傳統(tǒng)的遺傳學(xué)基因定位方法往往需要構(gòu)建漸滲系（recombinant inbred line，RIL）群體和復(fù)雜的遺傳圖譜［3-6］，但這種方式會(huì)耗費(fèi)大量的人力和時(shí)間成本。

分離群體分組分析方法（bulked segregant analysis，BSA）常用于鑒定與突變體表型緊密連鎖的基因區(qū)域或開發(fā)連鎖的分子標(biāo)記［7］。在傳統(tǒng)BSA方法中，得到突變體后，將突變體與親本雜交，在后代群體中根據(jù)突變表型對(duì)個(gè)體分組得到兩個(gè)分離群體，分別構(gòu)建兩個(gè)具有表型差異的基因池。在不同基因池中，與突變表型連鎖的分子標(biāo)記會(huì)表現(xiàn)出多態(tài)性，由此可通過比較分子標(biāo)記在不同基因池中的差異實(shí)現(xiàn)對(duì)突變基因的定位［8］。近年來，隨著第二代測(cè)序技術(shù)（next-generation sequencing，NGS）的逐漸成熟，衍生出眾多將BSA與NGS相結(jié)合的基因組定位策略BSA-seq，如MutMap［9］、QTLseq［10］、QTG-seq［11］和GradedPool-Seq（GPS）［12］等。這些方法大大縮短了基因定位的周期，同時(shí)降低了成本。

目前，BSA-seq定位策略已被應(yīng)用于挖掘黃瓜［13］、水稻［14］、小麥［15］、玉米［16］、花生［17］、大豆［18］等農(nóng)作物的優(yōu)良基因。程鳳等［19］利用MutMap 定位策略，以黃瓜中短果突變體為試驗(yàn)材料，在黃瓜1號(hào)染色體中發(fā)現(xiàn)影響果長(zhǎng)的基因CsaV3_1G044310。在小麥抗條銹病研究中，研究者們利用BSA-seq 定位到多個(gè)相關(guān)的數(shù)量性狀基因座（quantitative trait locus，QTL）位點(diǎn)［20-22］。

輻射誘變具有簡(jiǎn)單快捷、突變效率高、突變類型多樣等優(yōu)勢(shì)，常用于正向遺傳學(xué)研究。在已有的報(bào)道中，采用BSA-seq定位策略的研究中突變體材料大多來源于自然突變或化學(xué)誘變，對(duì)輻射誘變得到的突變體采用BSA-seq 方法定位突變基因的報(bào)道較少。as2-7D突變體是ASYMMETRIC LEAVES（AS2）基因的突變體［23］。AS2 基因編碼一個(gè)含有LATERAL ORGAN BOUNDARIES（LOB）結(jié)構(gòu)域的轉(zhuǎn)錄因子，該蛋白質(zhì)與ASYMMETRIC LEAVES 1（AS1）蛋白互作形成蛋白復(fù)合體，共同參與調(diào)節(jié)擬南芥葉片所有三個(gè)軸的正常形態(tài)［24］。as2-7D突變體的萼片背面（遠(yuǎn)離內(nèi)部花器官的一面）具有凸起，為明顯的表型特征，易于進(jìn)行表型抑制突變體的篩選。鑒于此，本研究以擬南芥突變體as2-7D為材料，對(duì)其進(jìn)行輻射誘變，從后代中篩選萼片背面較為平整光滑的抑制突變體。通過MutMap、QTL-seq 和GPS 三種BSA-seq 方法對(duì)其輻射誘變產(chǎn)生的突變體的基因進(jìn)行定位鑒定，探索BSA-seq 方法用于輻射誘變產(chǎn)生的突變體的基因定位的可能，旨在為加快輻射誘變基因的定位過程提供新思路。

1 材料與方法

1.1 試驗(yàn)材料

擬南芥野生生態(tài)型Columbia-0 （Col-0）和Landsbergerecta（Ler）由浙江大學(xué)原子核農(nóng)業(yè)科學(xué)研究所洪麗蘭課題組（本課題組）保存。擬南芥突變體as2-5D和SALK_127850 訂購(gòu)于擬南芥生物資源保藏中心（Arabidopsis Biological Resource Center，ABRC），as2-7D為本課題組前期篩選得到的材料，as2-7D和as2-5D具有相同的點(diǎn)突變，均在AS2 基因啟動(dòng)子上游1 484 bp 處發(fā)生了一個(gè)G 到A 的點(diǎn)突變，造成AS2 基因的異位表達(dá)［23］。as2-7D是Ler 背景，as2-5D是Col-0背景。因?yàn)閍s2-7D比as2-5D具有更顯著的萼片背面凸起的表型，故以as2-7D為材料進(jìn)行抑制突變體篩選。

取as2-7D自交系種子，在浙江省農(nóng)業(yè)科學(xué)院作物與核技術(shù)利用研究所用60Co-γ 射線誘變，輻射劑量率10 Gy·min-1，輻射劑量600 Gy，M1代植株自交一代，在M2代植株群體中篩選到若干萼片背面凸起表型顯著受到抑制的突變體，選取其中一個(gè)抑制突變體進(jìn)行下一步試驗(yàn)。

1.2 表型觀察

試驗(yàn)所用擬南芥的生長(zhǎng)溫度為22 ℃、光周期為16 h 光/8 h 暗，光照強(qiáng)度約為100 μmol·m-2·s-1。植物開花后一周，在DF295 體視鏡（Leica，德國(guó)）下觀察比較擬南芥花萼的表型差異并拍照。

1.3 遺傳學(xué)分析

將抑制突變體與Ler 和as2-7D分別進(jìn)行雜交，觀察F1代表型并分析突變基因的顯隱性；F1代自交得到F2，統(tǒng)計(jì)F2代植株群體的表型分離比，通過卡方測(cè)驗(yàn)分析抑制突變體萼片凸起受到抑制的表型是否由單基因所調(diào)控，從而確定突變基因的遺傳模式。

1.4 BSA-seq定位抑制基因

1.4.1 混池測(cè)序 取100 株親本as2-7D的幼嫩葉片構(gòu)建親本混池，為親本池。從抑制突變體和as2-7D雜交得到F2代群體中挑選as2-7D表型和抑制表型的植株各100株，取幼嫩葉片構(gòu)建子代混池，分別為as2-7D池和抑制突變體池。使用十六烷基三甲基溴化銨（cetyltrimethylammonium bromide，CTAB）提取液分別對(duì)三個(gè)混池樣本提取DNA，送往北京諾禾致源科技股份有限公司測(cè)序，DNA 樣本質(zhì)檢合格后構(gòu)建文庫，在NovaSeq 6000 測(cè)序平臺(tái)（Illumina，美國(guó)）進(jìn)行雙端重測(cè)序，測(cè)序深度為30×。使用fastp v0.23.1 軟件［25］對(duì)原始測(cè)序結(jié)果進(jìn)行過濾篩選，得到凈數(shù)據(jù)（clean data）。

1.4.2 測(cè)序結(jié)果分析 對(duì)測(cè)序結(jié)果的分析在服務(wù)器上進(jìn)行，服務(wù)器的操作系統(tǒng)為Ubuntu 20.04.4 LTS。從1001Genomes（https：//1001genomes. org/data/MPIPZ/MPIPZJiao2020/releases/current/strains/Ler/）上下載野生型Ler 的參考基因組序列以及基因注釋信息文件。使用fastqc v0.11.8 軟件查看測(cè)序質(zhì)量，質(zhì)量合格后分別使用MutMap、QTL-seq和GPS分析測(cè)序結(jié)果。

使用MutMap 方法分析測(cè)序結(jié)果：從github 上下載由Sugihara等［26］開發(fā)的MutMap pipeline（https：//github.com/YuSugihara/MutMap），使用該pipeline 默認(rèn)參數(shù)和實(shí)際混池樣本數(shù)，分析親本池和抑制突變體池的clean data，得到抑制基因的候選區(qū)間。

使用QTL-seq 方法分析測(cè)序結(jié)果：從github 上下載由Sugihara 等［26］開發(fā)的QTL-seq pipeline（https：//github.com/YuSugihara/QTL-seq），使用該pipeline 默認(rèn)參數(shù)和實(shí)際混池樣本數(shù)，分析親本池、as2-7D池和抑制突變體池的clean data，得到抑制基因的候選區(qū)間。

使用GPS方法分析測(cè)序結(jié)果：利用GPS方法對(duì)as2-7D池和抑制突變體池的clean data進(jìn)行分析。首先利用bwa v0.7.17軟件［27］將clean data與Ler的參考基因組比對(duì)定位。使用Picard軟件（https：//broadinstitute.github.io/picard）標(biāo)記重復(fù)序列。利用GATK v4.2.4.0軟件［28］實(shí)現(xiàn)對(duì)單核苷酸多態(tài)性（single nucleotide polymorphisms，SNP）和插入缺失（insertion-deletion mutations，InDel）變異位點(diǎn)的檢測(cè)。過濾掉低質(zhì)量和低測(cè)序深度的SNP和InDel 位點(diǎn)，使用Ridit（relative to an identified distribution unit）分析法計(jì)算-lg（P）值，得到與參考基因組序列具有顯著差異的突變位點(diǎn)及抑制基因的候選區(qū)間。

1.4.3 篩選可信的突變位點(diǎn) 在BSA-seq 的分析流程中得到基因序列變異信息（variant call format，vcf）文件，使用SnpEff v4.3t 軟件［29］對(duì)所有突變位點(diǎn)進(jìn)行注釋，結(jié)合三種BAS-seq 分析策略和Integrative Genomics Viewer （IGV）軟件［30］的基因組可視化結(jié)果，實(shí)現(xiàn)對(duì)抑制基因的精細(xì)定位。

2 結(jié)果與分析

2.1 as2-7D抑制突變體的獲得及表型特征

在as2-7D的輻射誘變?nèi)后w中篩選到一株萼片背面凸起表型顯著受到抑制的突變體。野生型Ler 萼片背面光滑平整（圖1-A），as2-7D萼片背面具有大量的凸起和褶皺（圖1-B），而抑制突變體的萼片背面較為平整，沒有凸起和褶皺，更接近于Ler野生型（圖1-C）。對(duì)該突變體的AS2基因測(cè)序，發(fā)現(xiàn)其仍然含有as2-7D的突變位點(diǎn)，且該基因區(qū)域不具有第二個(gè)突變位點(diǎn)，說明該突變體是由AS2基因之外的位點(diǎn)突變導(dǎo)致的對(duì)as2-7D萼片表型的抑制，選擇該抑制突變體進(jìn)一步研究。

圖1 as2-7D ssp1抑制突變體的表型Fig.1 The phenotypes of as2-7D ssp1 mutant

2.2 as2-7D抑制突變體的遺傳分析

將抑制突變體分別與as2-7D和Ler 雜交，觀察F1代的表型，抑制突變體 ×as2-7DF1代的萼片表型與as2-7D萼片表型一致（圖1-B、F），說明抑制突變體對(duì)as2-7D萼片表型的抑制是由隱性突變控制的。統(tǒng)計(jì)F1代自交得到的F2代植株群體的表型分離比，群體共43 個(gè)植株，其中類as2-7D32 株，類抑制突變體11 株，經(jīng)卡方測(cè)驗(yàn)分析得P值約為0.93，大于0.05，符合3∶1的分離比（表1），說明抑制突變體抑制as2-7D萼片表型的性狀是由一個(gè)單基因控制的，將該抑制基因命名為SUPPRESSORSofas2-7D SEPAL PHENOTYPE 1（SSP1），由此抑制突變體命名為as2-7D ssp1。

表1 as2-7D ssp1 × as2-7D F2代表型分離統(tǒng)計(jì)表Table 1 Phenotypic segregation of the as2-7D ssp1 × as2-7D F2 population

as2-7D ssp1× Ler F1代植株萼片背面仍然出現(xiàn)凸起，但凸起的數(shù)量較as2-7D有所減少（圖1-E），與as2-7D× Ler F1代具有相似的表型（圖1-D）。由于該表型凸起的數(shù)量介于Ler 和as2-7D之間，將該表型稱為中間型。as2-7D ssp1× Ler F2代群體中，類as2-7D∶中間型∶類Ler∶類as2-7D ssp1∶新表型的個(gè)體數(shù)量比符合3∶6∶3∶3∶1 的分離比（表2）。綜合上述試驗(yàn)數(shù)據(jù)推測(cè)，ssp1是一個(gè)單基因隱性突變，且不與AS2基因連鎖；as2-7D ssp1× Ler F2代群體中的新表型可能由SSP1單基因突變導(dǎo)致。

表2 as2-7D ssp1 × Ler F2代表型分離統(tǒng)計(jì)表Table 2 Phenotypic segregation of the as2-7D ssp1 × Ler F2 population

2.3 ssp1突變基因定位

對(duì)as2-7D親本池及as2-7D ssp1×as2-7DF2代群體的兩個(gè)表型池進(jìn)行全基因組測(cè)序，分別用MutMap、QTL-seq、GPS 方法分析測(cè)序結(jié)果，限定SSP1基因所在的位置，同時(shí)用基因組可視化軟件IGV 篩選候選基因。

圖2、3 分別為MutMap 和QTL-seq 的分析結(jié)果，藍(lán)色散點(diǎn)表示某些突變位點(diǎn)的SNP-index 值，綠色折線表示95%置信區(qū)間的閾值，黃色折線表示99%置信區(qū)間的閾值，紅色折線為2 Mb 內(nèi)的平均SNP-index。MutMap 法引入了SNP-index，表示混池中非參考基因組類型的SNP頻率。抑制突變體池中的某SNP標(biāo)記與抑制表型的關(guān)聯(lián)度越高，該標(biāo)記的SNP-index 越接近于1。查看MutMap 分析結(jié)果圖，在2 號(hào)染色體大概9 Mb 以及13～16 Mb 的位置范圍內(nèi)，平均SNP-index 已超過99%的置信閾值，但是在9～11 Mb 區(qū)間內(nèi)無可信的SNP 位點(diǎn)，未能繪制折線，使得后續(xù)2 Mb 范圍內(nèi)的平均SNP-index 值偏低。因此，根據(jù)MutMap 的結(jié)果，將SSP1基因限定在2號(hào)染色體9～16 Mb區(qū)間內(nèi)（圖2）。理論上，as2-7D ssp1×as2-7DF2代群體中兩個(gè)極端表型混池在SSP1基因所在區(qū)域的SNP-index 具有顯著差異，而其他區(qū)域差異較小。因此，可使用QTL-seq定位SSP1基因。QTL-seq 分別計(jì)算子代兩個(gè)極端表型混池的SNP-index，將兩個(gè)混池的SNP-index 做減法得到?SNP-index。?SNP-index越高的基因區(qū)域，與SSP1基因的關(guān)聯(lián)度也越高。從QTL-seq 分析圖同樣得出SSP1基因位于2號(hào)染色體9～16 Mb區(qū)間內(nèi)的結(jié)論（圖3）。GPS 分析結(jié)果顯示，2號(hào)染色體9～15 Mb區(qū)間內(nèi)富集了高-lg（P）值的散點(diǎn)（圖4），與MutMap 和QTL-seq 的分析結(jié)果基本一致。

圖2 MutMap分析定位SSP1基因的結(jié)果Fig.2 Results of MutMap analysis to locate the SSP1 gene

圖3 QTL-seq分析定位SSP1基因的結(jié)果Fig.3 Results of QTL-seq analysis to locate the SSP1 gene

圖4 GPS分析定位SSP1基因的結(jié)果Fig.4 Results of GPS analysis to locate the SSP1 gene

在2 號(hào)染色體的9～16 Mb 區(qū)間內(nèi)，共有391 個(gè)SNP和Indel 突變位點(diǎn)，根據(jù)質(zhì)量進(jìn)行初步篩選，得到139個(gè)可信的突變位點(diǎn)。以SNP-index 和ΔSNP-index 為指標(biāo)進(jìn)一步過濾候選位點(diǎn)，共得到2 個(gè)SNP-index 或ΔSNP-index高于99%置信度的候選位點(diǎn)。snpEff軟件的注釋信息表明這兩個(gè)突變位點(diǎn)位于基因的上游或下游，影響其附近基因功能的概率較低，排除這兩個(gè)突變位點(diǎn)為候選突變位點(diǎn)。

基因可視化軟件IGV 發(fā)現(xiàn)，抑制突變體混池在2 號(hào)染色體11 950 417～11 975 606 bp 范圍內(nèi)發(fā)生了一個(gè)25 189 bp的缺失（圖5）。為驗(yàn)證在as2-7D ssp1群體中存在該片段缺失，分別在缺失片段兩端設(shè)計(jì)鑒定引物，對(duì)6 株as2-7D ssp1植株的基因組DNA 進(jìn)行擴(kuò)增，結(jié)果顯示，as2-7D ssp1確實(shí)存在該大片段缺失。在該缺失片段中共包含AT2G28580、AT2G28590、AT2G28600、AT2G28605、AT2G28610、AT2G28620 這6 個(gè)基因。其中，AT2G28610編碼一個(gè)轉(zhuǎn)錄因子PRESSED FLOWER（PRS），在擬南芥?zhèn)壬鞴伲ㄈ巛嗥?、葉等）邊緣表達(dá)，具有調(diào)控細(xì)胞增殖的功能。根據(jù)前人的研究報(bào)道，PRS過表達(dá)會(huì)導(dǎo)致擬南芥萼片背面出現(xiàn)凸起的表型［31］，與as2-7D類似，因此將PRS作為SSP1基因的候選基因。

圖5 as2-7D ssp1突變體基因組中的大片段缺失Fig.5 Large deletion in the genome of the as2-7D ssp1

2.4 PRS作為抑制基因的驗(yàn)證

為驗(yàn)證PRS基因是否是SSP1基因，從Col-0 背景的PRS基因T-DNA 插入突變體SALK_127850，分離出純合的PRS基因突變體prs-3，構(gòu)建了as2-5D與prs-3的雙突變體as2-5D prs-3。觀察Col-0、as2-5D、prs-3和as2-5D prs-3 的萼片表型發(fā)現(xiàn)，as2-5D的萼片背面具有少量的凸起（圖6-B），prs-3 萼片背面平整光滑（圖6-C），與Col-0（圖6-A）相比無明顯差異，而as2-5D prs-3 的萼片背面沒有凸起，整體較為平整（圖6-D），說明PRS基因的缺失可以抑制as2-5D萼片背面凸起的表型，驗(yàn)證了PRS基因是SSP1基因。

圖6 PRS突變抑制as2-5D突變體的萼片表型Fig.6 PRS mutation suppresses the sepal phenotype of as2-5D mutant

3 討論

本研究使用60Co-γ 射線輻射誘變as2-7D，在輻射誘變后代中篩選得到抑制突變體as2-7Dssp1，成功地運(yùn)用BSA-seq 定位策略對(duì)抑制基因SSP1進(jìn)行了定位。BSA-seq 雖然無法精確定位大片段的插入或缺失突變，但可以在全基因組范圍內(nèi)篩選出于目標(biāo)基因高度連鎖的SNP 位點(diǎn)，從而將目標(biāo)基因的位置限制在較小的范圍內(nèi)。在得到目標(biāo)基因的候選區(qū)間后，可通過基因組注釋文件對(duì)該范圍內(nèi)的可信SNP 位點(diǎn)進(jìn)行篩選，并結(jié)合IGV軟件查看該區(qū)間內(nèi)是否存在大片段的插入或缺失，從而實(shí)現(xiàn)對(duì)目標(biāo)基因的精細(xì)定位。為使用BSA-seq 定位策略挖掘由輻射導(dǎo)致的突變基因提供了參考依據(jù)。

本研究通過MutMap、QTL-seq 和GPS 三種BSAseq 分析方法對(duì)輻射誘變突變體的基因進(jìn)行定位。MutMap 基于子代突變體池的SNP-index 對(duì)目標(biāo)性狀的基因區(qū)域進(jìn)行定位。在子代突變體群體的遺傳背景與親本一致或高度相似時(shí)，MutMap具有良好的定位效果。本研究中的抑制突變體由as2-7D輻射誘變得到，二者的雜交F2代群體具有相似的遺傳背景，因此可用MutMap定位抑制基因。QTL-seq的技術(shù)思路與MutMap相類似。為減小背景噪音的影響，對(duì)兩個(gè)子代混池分別求SNP-index，相減后得?SNP-index，根據(jù)?SNPindex 對(duì)抑制基因進(jìn)行定位。本研究中MutMap 和QTL-seq 均將抑制基因定位在2 號(hào)染色體的9～16 Mb區(qū)間內(nèi)。GPS 方法沒有使用SNP-index 算法，而是使用Ridit 分析法來篩選與目標(biāo)基因關(guān)聯(lián)度較高的SNP位點(diǎn)，通過窗口smooth 來減小背景噪音的影響。本研究中，三種BSA-seq 分析方法均成功定位到了抑制突變體的變異區(qū)間，且取得了相似的定位效果。因此，MutMap、QTL-seq 和GPS 均可應(yīng)用于輻射誘變突變體的基因定位。

4 結(jié)論

本研究在as2-7D的輻射誘變突變體后代中篩選到抑制突變體as2-7D ssp1，使用MutMap、QTL-seq 以及GPS 三種BSA-seq 定位策略對(duì)引起突變性狀的SSP1基因進(jìn)行定位，三種方法均成功定位到了抑制突變體的變異區(qū)間；并通過基因組注釋信息和遺傳學(xué)試驗(yàn)分析，確認(rèn)變異區(qū)間內(nèi)的PRS基因的缺失導(dǎo)致了as2-7D ssp1的突變表型。結(jié)果表明BSA-seq定位策略可用于定位輻射誘變引起的大片段缺失突變。