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

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

（1浙江大學原子核農業科學研究所/農業農村部和浙江省核農學重點實驗室，浙江 杭州 310058；2浙江大學生命科學學院植物生物學研究所，浙江 杭州 310058）

從人工誘導或自然發生的突變體中挖掘出突變信息既是分子遺傳研究的重要內容，也是培育高產優質新品種的重要途徑［1-2］。傳統的遺傳學基因定位方法往往需要構建漸滲系（recombinant inbred line，RIL）群體和復雜的遺傳圖譜［3-6］，但這種方式會耗費大量的人力和時間成本。

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

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

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

1 材料與方法

1.1 試驗材料

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

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

1.2 表型觀察

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

1.3 遺傳學分析

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

1.4 BSA-seq定位抑制基因

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

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

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

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

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

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

2 結果與分析

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

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

圖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），說明抑制突變體對as2-7D萼片表型的抑制是由隱性突變控制的。統計F1代自交得到的F2代植株群體的表型分離比，群體共43 個植株，其中類as2-7D32 株，類抑制突變體11 株，經卡方測驗分析得P值約為0.93，大于0.05，符合3∶1的分離比（表1），說明抑制突變體抑制as2-7D萼片表型的性狀是由一個單基因控制的，將該抑制基因命名為SUPPRESSORSofas2-7D SEPAL PHENOTYPE 1（SSP1），由此抑制突變體命名為as2-7D ssp1。

表1 as2-7D ssp1 × as2-7D F2代表型分離統計表Table 1 Phenotypic segregation of the as2-7D ssp1 × as2-7D F2 population

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

表2 as2-7D ssp1 × Ler F2代表型分離統計表Table 2 Phenotypic segregation of the as2-7D ssp1 × Ler F2 population

2.3 ssp1突變基因定位

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

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

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

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

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

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

基因可視化軟件IGV 發現，抑制突變體混池在2 號染色體11 950 417～11 975 606 bp 范圍內發生了一個25 189 bp的缺失（圖5）。為驗證在as2-7D ssp1群體中存在該片段缺失，分別在缺失片段兩端設計鑒定引物，對6 株as2-7D ssp1植株的基因組DNA 進行擴增，結果顯示，as2-7D ssp1確實存在該大片段缺失。在該缺失片段中共包含AT2G28580、AT2G28590、AT2G28600、AT2G28605、AT2G28610、AT2G28620 這6 個基因。其中，AT2G28610編碼一個轉錄因子PRESSED FLOWER（PRS），在擬南芥側生器官（如萼片、葉等）邊緣表達，具有調控細胞增殖的功能。根據前人的研究報道，PRS過表達會導致擬南芥萼片背面出現凸起的表型［31］，與as2-7D類似，因此將PRS作為SSP1基因的候選基因。

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

2.4 PRS作為抑制基因的驗證

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

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

3 討論

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

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

4 結論

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