基于BSA-seq結合連鎖分析發掘大豆莢粒性狀QTLs及候選基因

孫亞倩， 陳士亮， 褚佳豪， 李喜煥*， 張彩英*

（1.河北農業大學，華北作物改良與調控國家重點實驗室，華北作物種質資源研究與利用教育部重點實驗室，河北 保定 071001；2.承德市農林科學院，河北 承德 067000）

大豆是重要糧油作物，提供了全球50%以上的食用油和25%以上的植物蛋白[1]。然而，大豆單產水平和單位面積經濟效益較其他農作物而言相對較低，亟需進一步提高以滿足日益增長的需求。豆莢和籽粒是菜用和粒用大豆重要的產量相關性狀，直接影響單位面積產量和經濟效益[2]。莢粒性狀還是菜用和粒用大豆最重要的外觀品質性狀，對其商品價值具有重要影響。因此，研究莢粒性狀對提高大豆產量、品質、商品價值和經濟效益等具有重要意義。

大豆莢粒性狀包括莢長、莢寬、莢重、粒長、粒寬和粒重等，這些性狀屬于數量性狀，且已獲得少數控制其遺傳位點及候選基因。Niu等[3]獲得59個粒長、粒寬和粒厚等性狀QTLs，位于20條染色體；Xie 等[4]通過精細定位6 號染色體籽粒大小QTLs獲得8個候選基因；Dargahi等[5]獲得16個粒長、粒寬和粒厚QTLs，位于9 條染色體，表型貢獻率7.4%~18.4%；Yang 等[6]檢測到60 個籽粒性狀QTLs，位于除5 和19 號以外的18 條染色體；Teng等[7]研究發現，5 個粒長QTLs 位于7、12、13 和17號染色體，3 個粒厚QTLs 位于9、12 和13 號染色體，5 個粒寬QTLs 位于12、13 和17 號染色體；Cui等[8]獲得30 個籽粒性狀QTLs，位于12 條染色體，包括1 個多環境下的一因多效QTL；Hina 等[9]認為，6、10、13和20號染色體存在籽粒性狀QTLs 熱點區；Kumawat 等[10]獲得53 個籽粒大小QTLs，位于除3、14和18號以外的17條染色體。

目前，隨著測序技術的發展，混合群體分離測序（bulked segregant analysis sequencing，BSA-Seq）因無需構建分子遺傳連鎖圖譜，在重要性狀遺傳位點發掘中越來越受到重視[11?12]。張之昊等[13]采用BSA-seq 技術挖掘大豆的多小葉基因，經SNPindex 算法獲得11號染色體的3個關聯區域，總長度3.42 Mb，包含701 個基因；曾維英等[14]采用BSA-seq 技術挖掘大豆抗豆卷葉螟基因，經SNPindex 算法獲得329 個基因，位于7、16 和18 號染色體；刁衛楠等[15]通過BSA-seq 技術將控制西瓜黃色果肉的主效QTLs 定位于6 號染色體，并結合基因注釋和表達分析獲得5 個候選基因；Zhao等[16]利用BSA-seq 方法發掘到水稻3 號染色體控制粒長的基因熱點區，并通過基因編輯技術證實了粒長基因；Zhang 等[17]利用BSA-seq 方法獲得3 個水稻株高QTLs、2 個穗長QTLs 和4 個抽穗期QTLs。

綜上，盡管目前已有關于大豆莢粒性狀遺傳位點發掘與相關基因研究，但多基于二代分子標記連鎖定位方法獲得，因圖譜密度和飽和度較低，所得連鎖標記及基因數目少，且很難在育種實踐中應用。鑒于此，本研究基于大豆重組自交系群體（recombinant inbred line, RIL）在多種環境條件下的莢粒性狀（莢長、莢寬、莢重、粒長、粒寬和粒重），篩選極端家系并構建2個混池，通過BSA-seq技術發掘控制莢粒性狀的遺傳位點，結合前期通過連鎖定位方法獲得的該群體QTLs[18]，進一步發掘在多環境條件下、2種方法同時能檢測到的一因多效遺傳位點；同時結合RIL群體親本莢粒不同發育時期轉錄組數據篩選候選基因，為開展大豆莢粒性狀分子遺傳改良和遺傳基礎解析提供支撐。

1 材料與方法

1.1 試驗材料

供試大豆重組自交系群體（C813×KN7，F6:8，193個家系）[18]由河北農業大學大豆遺傳育種課題組提供，其親本C813為大莢大粒菜用大豆優異種質（鮮莢長 5.0 cm、鮮莢寬1.3 cm、鮮莢重2.2 g、鮮粒長1.4 cm、鮮粒寬0.8 cm、鮮粒重0.6 g），KN7 為小莢小粒育成品種（鮮莢長4.2 cm、鮮莢寬1.1 cm、鮮莢重1.4 g、鮮粒長1.2 cm、鮮粒寬0.8 cm、鮮粒重0.4 g）。

1.2 試驗方法

1.2.1大豆重組自交系群體種植 供試大豆重組自交系群體分別于2020、2021 年種植在4 種環境條件下，即2021 年保定（E1）、2020 年保定（E2）、2020 年石家莊鹿泉（E3）、2020 年廊坊文安（E4），田間試驗采用完全隨機區組設計，行長3 m，行距0.5 m，株距10 cm，3次重復，田間種植和管理參考Chen等[18]方法。

1.2.2大豆重組自交系群體莢粒性狀鑒定 待大豆重組自交系群體生長至R6~R7 期，取植株中部生長一致的標準二粒莢10 個，測定其莢長（pod length, PL）、莢寬（pod width, PW）和莢重（pod weight, PWT），重復3 次[18]；隨后剝去莢皮，將所有籽粒置于SC-G種子自動分析系統（杭州萬深檢測科技有限公司）樣品板上，分析粒長（seed length,SL）、粒寬（seed width, SW）和粒重（seed weight,SWT）[2,18-20]。

1.2.3重組自交系群體莢粒性狀極端家系篩選根據供試大豆重組自交系群體各家系在不同環境條件下莢粒性狀的綜合表現，分別選擇莢長>4.5 cm且莢寬>1.2 cm的20個大莢大粒極端家系和莢長<4.5 cm 且莢寬<1.2 cm 的20個小莢小粒極端家系，用于BSA-seq發掘控制大豆莢粒性狀QTLs。

1.2.4大豆莢粒性狀混池BSA-seq 分別提取上述入選的家系以及RIL群體親本（C813和KN7）基因組DNA，并將其按照性狀分類，等量混合成2個混池——大莢大粒池（large pool，LP）和小莢小粒池（small pool，SP），分別含有20 個家系。將LP、SP、C813 和KN7 的DNA 送廣州基迪奧生物科技有限公司完成BSA-seq 分析，測序深度為30×，參考基因組為Williams82 v2.0，測序質量通過Q20（高通量測序中，錯誤率低于1%的堿基百分比）和Q30（高通量測序中，錯誤率低于0.1%的堿基百分比）等進行評估。

另外，為保證結果可靠性，按照以下標準篩選SNP 和Indel 標記：①入選標記在親本間存在差異，且符合RIL 群體分離方式；②RIL 群體親本測序深度大于閾值；③標記在2 個混池中均不缺失；④每個混池測序深度需大于10×且小于500×；⑤至少1 個混池的SNP 指數（SNP index，SNPI）在0.3~0.7之間。

1.2.5發掘莢粒性狀遺傳位點的算法 本研究采用3 種算法進行莢粒性狀遺傳位點發掘，包括SNP index 算法、歐氏距離（Euclidean distance,ED）算法和G-statistic 算法[13,21]；置信水平分別設為0.95和0.99，3種算法的計算公式如下。

①SNP index算法。

式中，SNPI為SNP 指數，ρ 為位點在混池中的深度。

②ED算法。

式中，mut與wt為隱性和顯性混池，A、C、G、T為各等位標記位點的覆蓋深度。

③G-statistic算法。

式中，O為等位深度的觀測值，E為等位深度的期望值，ln為自然對數，i為等位基因數目。

1.2.6BSA-seq 一致性區間候選基因篩選 首先，比較供試大豆RIL群體BSA-seq不同算法關聯結果，篩選一致性關聯區間；然后，在一致性區間內尋找候選基因，并結合基因注釋、基因DNA 水平SNP 等位變異以及RNA 水平表達分析，篩選大豆莢粒性狀相關基因。候選基因查找及其基因注釋參考大豆測序基因組Williams82（Glyma2.0）（https://www.soybase.org/），基因的RNA 水平表達分析采用供試RIL群體親本C813和KN7的豆莢5個不同發育時期（Pod-T1~Pod-T5）和籽粒3 個不同發育時期（Seed-T1~Seed-T3）的基因表達轉錄組數據[2]。

1.2.7表型數據分析 利用Microsoft Excel 2017

統計軟件的成組數據t檢驗方法，對供試RIL群體極端家系在單一環境條件下以及多種環境條件下的6 個莢粒相關性狀（莢長、莢寬、莢重、粒長、粒寬和粒重）進行顯著性檢驗。

2 結果與分析

2.1 大豆重組自交系群體BSA-seq 極端家系的遴選

通過分析供試大豆重組自交系群體2 年4 種環境條件下的莢長、莢寬、莢重、粒長、粒寬和粒重，篩選出40 份極端家系構成2 個混池——LP 和SP。分析2 個混池不同環境條件下的莢粒性狀（表1）發現，2 個混池在單一環境條件下的莢長、莢寬、莢重、粒長、粒寬和粒重均存在極顯著差異，并且在多環境條件下的莢粒性狀亦存在顯著或極顯著差異，為進一步利用BSA-seq 方法發掘控制其遺傳位點奠定了重要基礎。

表1 大豆RIL群體BSA-seq極端家系不同環境莢粒性狀遺傳差異Table 1 Genetic variation of extreme RIL lines on pod and seed traits for BSA-seq under different environments

2.2 BSA-seq測序質量及SNP與Indel分析

通過分析2 個混池及其親本材料的BSA-seq測序質量（表2）發現，過濾后的有效數據量為30 Gb，Q20 在97.2%以上，Q30 在92.5%以上，G和C 堿基含量為36.32%~36.58%，過濾后的總讀段為2 Gb；與參考基因組的平均比對效率在98%以上，且基因組覆蓋度深，4 個樣品的10×基因組覆蓋率均在94%以上，20×基因組平均覆蓋率則在80%以上，說明4 個樣品的BSA-seq 測序數據充足且質量高，可進行后續BSA-seq 關聯位點分析。

表2 兩個混池及其親本BSA-seq測序質量分析Table 2 Quality analysis of BSA-seq of two pools and soybean RIL parents

進一步分析BSA-seq 測序SNP（表3）與Indel（表4），共獲得162萬~241萬個SNP，其中位于外顯子區的有6.2 萬~9.1 萬個，基因上游的有9.6 萬~14.4萬個，引起非同義突變的有3.5萬~5.1萬個，引起提前終止的有726~1 089個，引起終止密碼子丟失的有144~223個；同時獲得30萬~43萬個Indel，其中位于外顯子區的有3 823~5 338個，基因上游的有3.0 萬~4.4 萬個，引起移碼缺失的有1 128~1 534個，移碼插入的有987~1 314個，引起提前終止的有69~95 個，引起終止密碼子丟失的有10~18個。

表3 兩個混池及其親本BSA-seq測序SNP數量與質量分析Table 3 SNP number and quality analysis of BSA-seq of two pools and soybean RIL parents

表4 兩個混池及其親本BSA-seq測序Indel數量與質量分析Table 4 Indel number and quality analysis of BSA-seq of two pools and RIL parents

2.3 大豆莢粒性狀BSA-seq遺傳位點發掘

通過分析篩選后的SNP 和Indel標記發現，符合入選標準的總標記數量為763 727個，其中包括SNP 標記658 186 個、Indel 標記105 541 個，并以17、9和3號染色體入選標記數量最多。在此基礎上，利用上述SNP 和Indel 標記分別通過SNPindex、ED 和G-statistic 算法發掘莢粒性狀關聯位點。

2.3.1SNP-index 算法發掘遺傳位點分析 采用SNP-index 算法分析與莢粒性狀關聯位點，結果（圖1A）表明，在置信度為0.95 時，獲得27 個關聯區域，位于15 條染色體（除5、8、10、11 和18 號外），其中7 號染色體包括7 個區域，16 和20 號染色體各包括3 個區域。同時發現，在置信度為0.99 時，獲得5 個關聯區域，分別位于5 條染色體（1、4、7、12和19號）。

圖1 大豆莢粒性狀不同算法關聯區域Fig. 1 Associated regions of soybean pod and seed traits via different methods

2.3.2ED算法發掘遺傳位點分析 采用ED方法分析與莢粒性狀關聯的位點，結果（圖1B）表明，在置信度0.95 時，獲得41 個關聯區域，分別位于17 條染色體（除5、8 和11 號外），其中7 和20 號染色體各包括5 個區域，6、10、13、14、16 和20 號染色體均包括3個區域；在置信度0.99時，獲得15個關聯區域，分別位于9 條染色體（1、3、4、7、12、13、15、19 和20 號），其中7 和20 號染色體包括6 個區域，4號染色體包括2個區域。

2.3.3G-statistic 算法發掘遺傳位點分析 采用G-statistic 方法分析與莢粒性狀關聯的位點，結果（圖1C）表明，在置信度0.95 時，獲得51 個關聯區域，位于18條染色體（除5和18號外），其中7號染色體包括10 個區域，20 號染色體包括8 個區域；在置信度0.99時，獲得16個關聯區域，分別位于9條染色體（1、3、6、7、12、15、17、19 和20 號），其中7 號染色體包括6 個區域，3 和20 號染色體包括2個區域。

2.3.43 種算法共關聯遺傳位點分析 綜合分析上述3 種算法獲得的關聯區域（表5），在置信度0.95 時，獲得27 個一致性關聯區域，位于15 條染色體（除5、8、10、11和18號外）；在置信度0.99時，獲得4 個一致性關聯區域，分別位于4 條染色體（1、7、12 和19 號），其中15 個基因位于1 號染色體，57 個基因位于7 號染色體，82 個基因位于12號染色體，182個基因位于19號染色體。

表5 大豆莢粒性狀不同算法一致性關聯位點Table 5 Genetic loci of soybean pod and seed traits via different methods

2.4 BSA-seq與連鎖定位一致性遺傳位點發掘

綜合分析本研究采用BSA-seq 方法獲得的遺傳位點以及課題組前期采用連鎖定位方法（SoySNP6K 分析）獲得的遺傳位點[18]發現，在7 號染色體的10.45~11.59 Mb 和12.96~13.54 Mb 區段，19 號染色體的47.63~48.50 Mb 和49.61~49.91 Mb 區段以及20 號染色體的35.83~36.24 Mb 區段存在一致性。并且發現，20 號染色體SNP 標記ss715637629（物理位置為35 836 361），19 號染色體SNP 標記ss715635705（物理位置為47 636 564）、ss715635755（物理位置為48 040 350）、ss715635827（ 物理位置為 48 375 196）、ss715635844（ 物理位置為 48 501 612）、ss715635790（ 物理位置為 48 196 242）、ss715635952（物理位置為49 618 939）和ss715635964（物理位置為49 699 850）在連鎖定位和BSA-seq分析中同時被檢測到，一方面證實了本研究結果的可靠性，另一方面為菜用和粒用大豆莢粒性狀分子遺傳改良提供了新的信息。

2.5 大豆莢粒性狀候選基因篩選

為進一步篩選大豆莢粒性狀候選基因，本研究在上述3 種算法0.99 置信水平共定位的染色體區段尋找到336 個莢粒性狀候選基因的SNP 等位變異（表6），有40 個基因在DNA 水平存在外顯子非同義突變，暗示其可能與大豆莢粒性狀有關。

表6 BSA-seq一致性區間存在外顯子非同義突變的候選基因Table 6 Candidate genes in the consistent regions with non-synonymous mutations in gene exons

在此基礎上，利用供試RIL 群體親本C813 和KN7的豆莢和籽粒不同發育時期轉錄組數據分析上述40 個基因，結果（圖2、圖3）發現，有15 個基因在豆莢和籽粒中表達量較高，其中包括Glyma.07G108700、Glyma.19G220000、Glyma.19G222200、Glyma.19G222300和Glyma.19G222400等；候選基因Glyma. 19G222200、Glyma. 19G222300和Glyma.19G224400隨豆莢發育進程表達量呈現增加趨勢。由此可見，上述15個候選基因在DNA水平存在外顯子非同義突變，而且RNA 水平在豆莢中表達量較高，故認為其與大豆莢粒性狀形成具有密切關系，尤其是隨著豆莢發育進程表達量呈現增加趨勢的候選基因Glyma.19G222200、Glyma.19G222300和Glyma.19G224400，可用于基因克隆以進一步證實其在大豆莢粒發育過程中的生物學功能。

圖2 大豆莢粒性狀候選基因在C813中的表達分析Fig. 2 Expressions of candidate genes related to soybean pod and seed traits in C813

圖3 大豆莢粒性狀候選基因在KN7中的表達分析Fig. 3 Expressions of candidate genes for soybean pod and seed traits in KN7

3 討論

3.1 大豆莢粒性狀一致性遺傳位點的發掘與利用

BSA策略可用于存在遺傳差異的多種群體，且可與傳統定位或關聯方法結合使用[21]。郭微[22]利用BSA-seq結合連鎖定位方法將水稻耐鹽堿QTLs縮小至2 號染色體的2 個區段；李君[23]利用BSAseq結合連鎖定位方法尋找玉米圓斑病遺傳位點；Vogel等[24]結合BSA-seq與連鎖定位方法獲得南瓜疫霉根腐病和冠腐病QTLs。本研究利用BSA-seq對大豆RIL 群體極端家系構建的2 個混池進行分析，經綜合3 種算法關聯分析，在0.95 置信水平將控制莢粒性狀遺傳位點定位在15 條染色體；同時，本課題組前期曾利用該RIL 群體，通過連鎖定位方法將控制莢粒性狀遺傳位點定位在7 條染色體，表型貢獻率5.28%~17.30%[18]。進一步比較BSA-seq 結果與連鎖定位結果發現，二者在7、19和20 號染色體存在共定位區間，且19 和20 號染色體的多個SNP 標記被2 種策略同時檢測到，充分證實了本研究結果可靠性和真實性。

關于20 號染色體控制大豆莢粒性狀遺傳位點，Hina 等[9]曾在該染色體的36.18~38.30 Mb 和35.92~38.14 Mb 區段檢測到控制粒長與長寬比QTLs；Zhao 等[25]曾在該染色體35.23 Mb 附近檢測到控制粒長QTLs；本課題組前期曾在該染色體35.84~36.23 Mb 區段檢測到多環境、多性狀共定位一因多效QTL（ss715637668~ss715637629），可同時控制大豆莢重、粒重、粒長和籽粒面積，解釋8.92%~17.30%表型變異，而且在相應性狀BLUP值的連鎖定位中同時也被檢測到，解釋12.07%~16.81%表型變異[18]。本研究通過BSA-seq 技術，將控制大豆莢粒性狀遺傳位點定位到20 號染色體的34.99~38.10 Mb 區段，與上述連鎖定位存在一致性區間。并且，為進一步利用上述一致性SNP標記，本課題組將該標記轉換為PCR標記，經擴增供試RIL 群體親本C813 和KN7，結果發現，能夠按照預先設計擴增出目的條帶（C813 可擴增500 bp 左右目的條帶，而KN7 中沒有該目的條帶），繼續擴增RIL 群體2 個混池（LP 和SP）的40個家系，發現LP 混池中有15 個家系可擴增出目的條帶，而SP 混池中僅有4 個家系可擴增出目的條帶，進一步擴增RIL群體其他家系將其劃分為2種基因型，且與預期結果基本一致，既證實了本研究結果的真實性，也為大豆莢粒性狀分子遺傳改良提供了技術支撐。

3.2 控制大豆莢粒性狀候選基因分析

關于大豆莢粒性狀候選基因，目前研究報道較少。Wang 等[1]獲得了可同時影響籽粒大小、脂肪和蛋白質含量的功能基因GmSWEET10a和GmSWEET10b，超表達2 個基因均可顯著增加百粒重和脂肪含量；而沉默基因后則可顯著降低百粒重和脂肪含量。最近，王明曉等[26]將控制豆科植物種子大小的基因劃分為3 類，即轉錄因子類（如WRKY、AP2/ERF 和TIFY 家族等）、植物激素類（生長素、油菜素內酯、赤霉素和細胞分裂素等）和其他因子類（細胞色素P450 家族、泛素介導的蛋白水解、細胞壁轉化酶抑制因子、肌醇單磷酸酶、ABC 轉運蛋白等），但同時也指出，豆科植物種子大小非常復雜，目前有關其控制基因及其調控網絡仍很欠缺。鑒于此，本研究在前述多年、多點鑒定RIL 群體6 個莢粒性狀，并分別通過BSAseq 和連鎖定位方法發掘遺傳位點的基礎上，進一步尋找控制莢粒性狀候選基因，結果發現，在一致性關聯區間內有336 個基因，并有40 個基因的外顯子存在非同義突變，其中包括轉錄因子基因Glyma. 19G221700、Glyma. 19G222200、Glyma.19G224700和Glyma.19G236900以及其他因子類基因Glyma.19G219500、Glyma.19G224400和Glyma.19G226000等。

在上述候選基因中，Glyma.19G222200編碼蛋白屬于MYB 轉錄因子，該基因在供試大豆RIL群體親本和混池間存在2 個非同義突變，且基因表達量隨大豆莢粒發育而呈現上升趨勢。擬南芥MYB 轉錄因子基因KUODA1參與植物葉片細胞的擴張，該基因突變后，可使轉基因植株的葉片顯著變小，而恢復植株的葉片與野生型基本相同，且基因超表達植株的葉片面積顯著增大[27]；擬南芥MYB 轉錄因子基因DRMY1參與植物細胞的擴張過程，該基因突變可引起一系列參與植物細胞壁生物合成與重塑的基因表達量發生改變，并使轉基因植株的角果顯著變短[28]。因此，本研究推測MYB 轉錄因子基因Glyma.19G222200可能通過影響其下游一系列與大豆莢粒生長發育相關基因的表達來控制莢粒發育，因而具有重要研究價值。

另外，植物體內的硝酸鹽轉運蛋白不僅具有運輸亞硝酸鹽的重要功能，而且可運輸生長素、赤霉素、脫落酸和茉莉酸等多種物質，因而對植物正常生長發育至關重要。魏一鳴[29]通過研究玉米EMS 突變體發現，硝酸鹽轉運蛋白基因ZmNPF7.9發生單堿基的非同義突變（G/A 突變），可導致玉米籽粒變小，頂部凹陷且粒重下降；進一步分析發現，該突變體植株參與糖酵解、碳代謝及氨基酸合成的關鍵基因表達量降低，推測其通過影響籽粒中的營養物質運輸和代謝進而導致玉米籽粒變小、粒重下降。在本研究中，候選基因Glyma.19G224400編碼蛋白屬于硝酸鹽轉運蛋白家族，該基因在供試大豆RIL 群體的2 個親本間以及混池間存在3 個非同義突變，且基因表達量隨大豆莢粒發育進程呈現上升趨勢，因而推測其可能通過參與大豆莢粒中營養物質以及激素類物質運輸進而影響莢粒發育。綜上，本研究認為，上述候選基因在DNA 水平存在外顯子SNP 非同義突變，而且隨大豆莢粒生長發育呈現上調表達，說明其與莢粒生長發育密切相關，為進一步研究基因生物學功能奠定了重要基礎。