BSA-Seq 結合RNA-Seq 技術挖掘大豆葉片提前黃化衰老基因

李世寬 洪慧龍 付佳祺 谷勇哲 孫如建 邱麗娟,*

1 哈爾濱師范大學生命科學與技術學院, 黑龍江哈爾濱 150025; 2 中國農業科學院作物科學研究所, 北京 100081; 3 呼倫貝爾市農業科學研究所, 內蒙古扎蘭屯 162650

葉片光合作用是植物生長發育的重要過程, 是農作物產量構成因素之一。葉片變黃是植物開始衰老的顯著特征并導致作物減產[1-5]。葉色黃化突變體相關研究始于1933 年, 陸地棉黃葉突變體[6]導致棉花減產, 甚至使棉花死亡。1948 年, Huang 等[7]利用X 射線誘變獲得小球藻失綠突變體wcm60, 驗證了原卟啉是葉綠素合成的前體, 葉色變異開始逐漸得到重視。

大豆葉片黃化衰老是物質從葉片向發育種子移動的結果[8]。大豆粒重絕大部分來源于開花后的光合產物[9]。在大豆籽粒形成期, 供應中心是本葉腋的籽粒[10]。鼓粒期葉片的光合作用對籽粒生長發育影響最大、其次是莢、再次是葉柄[11], 因此葉片提前黃化衰老會影響大豆產量。1997 年, Lin 等[12]利用Pride B216 × A15 和Anoka × A7 群體定位了缺鐵性黃化的7 個QTL, 揭示了缺鐵性黃化的多基因控制現象。

Palmer 等發現了非致死型黃化突變體T323, 及其等位的致死型黃化突變體T362H與T225H。用SSR 分子標記對Minsoy ×T225H群體作圖, 定位了致死型黃化表型的基因座y18_1與y18_2[13-14]。利用‘波高’ב南農94-156’的重組自交系將與大豆葉綠素累積量相關聯的4 個QTL 分別定位于D1a、H、M、G 和F 連鎖群[15]。近些年來也不斷有報道利用突變體和圖位克隆等方法, 將葉片黃化基因定位在1 號、2 號、5 號、10 號、11 號、13 號、15 號、18 號和20 號染色體[16-25], 廣泛存在于多條染色體上, 這說明導致葉片黃化的因素有很多。

基因挖掘技術不斷創新, 突變體是進行正向遺傳學研究的優異材料, 對隨機誘變產生的突變體進行研究, 不僅能夠挖掘新的基因, 也可以揭示已知基因的新功能[26-29]。混合分組分析法(Bulked Segregant analysis sequencing, BSA-Seq)在挖掘新基因研究中的應用日益廣泛[29-33], 該方法在不構建遺傳圖譜的情況下也能高效快速地挖掘目的基因。處理海量的RNA 測序(RNA sequencing, RNA-Seq)數據,挖掘出有用的信息逐漸成為了研究人員新目標。RNA-Seq 不僅能檢測基因表達量, 還可以檢測基因的等位變異, 是降低基因組復雜性的一種經濟有效的方法。雖然基于RNA-Seq 的SNP/InDel 挖掘只能檢測來自轉錄區域的變異信息, 但是大多數的來自轉錄區域的SNP/InDel 才是直接影響功能的變異,并且基于RNA-Seq 挖掘特定區域(如轉錄區域)的SNP/InDel 能提供非常高深度的變異信息, 還能檢測到高表達基因中低頻的變異[34-35]。競爭性等位基因特異性PCR (kompetitive allelespecific PCR, KASP)技術是基于PCR 的高通量SNP 分型技術, 可基于特定群體的變異信息, 在目標區域開發特異標記, 將圖位克隆與KASP 標記的連鎖分析相結合的策略,在作物基因挖掘研究中已得到應用[36-38]。

關于全生育期或營養生長期的葉色突變體報道較多, 而在大豆鼓粒后期出現葉片提前黃化衰老性狀尚未見報道。本研究在栽培大豆中品661 突變體庫中鑒定出一個大豆鼓粒后期非致死型葉片提前10～20 d 黃化突變體ly。利用BSA-Seq 方法初定位并圖位克隆, 結合RNA-Seq 技術確定了4 個SNP 和8 個差異表達候選基因, 為探索大豆鼓粒后期葉片提前黃化衰老圖位克隆和衰老分子機制研究提供新的思路。

1 材料與方法

1.1 試驗材料

在栽培大豆中品661 的EMS 誘變突變體庫中鑒定出一個大豆鼓粒后期葉片提前黃化突變體ly。利用ly與大豆鼓粒后期葉片正常色的野生型ofc配制組合得到F0, 在內蒙古扎蘭屯市農牧科學研究所溫室繁種, 獲得F1單株。在北京順義種植F2獲得978個單株, 2020 年種植 F2:3, 根據 F2:3表型取樣做BSA-Seq, F2:4取樣測RNA-Seq, 2022 年每株F2分出50 粒種子為株行, 海南調查F2:3, 鼓粒后期分別觀察F2:3株行表型并取樣。

1.2 表型分析

1.2.1 遺傳分析 觀察F2:3鼓粒后期每行的葉色, 將其分為黃色純合行、綠色純合行和分離雜合行, 其中低于30 株的株行表型數據不計入。用適合性測驗(χ2測驗)檢測F2:3黃純合行、綠純合行和分離行的行數比。假定無效假設為: F2:3群體的行數比符合1∶2∶1, 若χ2<5.99, 則P>0.05 接受無效假設, 若χ2>5.99, 則P<0.05 否定無效假設。計算公式如下:

式中,O為實際行數;E為理論行數。

1.2.2 生理指標測定 在F2:3鼓粒后期分離行中,分別選擇野生型和突變型表型的材料, 取鼓粒充分的3 粒莢, 每組取5 個生物重復。為排除時間和種植位置的影響, 在不同時間和株行分別進行8 次重復。稱莢果重和種子鮮重, 種子風干后稱干重, 用(種子鮮重-種子干重)/莢果重來表示水分, 種子干重/莢果重來表示干重, 其中莢果重用來消除莢的生長差異, 對得到的數據進行配對t檢驗, 若P<0.05則代表水分或干重有顯著差異; 在表型出現前大約1 周開始用SPAD-502 測定葉片葉綠素含量的相對值即SPAD 值, 直到葉片凋落, 每個性狀測3 株, 測3個不同位置葉片, 每個葉片測3 次SPAD 值取平均后計數, 對不同表型得到的數據去掉最大值和最小值后取平均, 并對不同時間的SPAD 值進行t檢驗若P<0.05 則代表此時的SPAD 值有顯著差異。

1.2.3 農藝性狀數據分析 在海南三亞冬季短日照條件下大田種植F2:3, 由于野生型親本ofc(矮皺葉, 圖1-A)和突變型親本ly(高展葉, 圖1-C)有明顯不同, 因此調查時分為矮皺葉和高展葉2 組, 且選擇生長均勻且長勢相近的純合行, 挑選F2:3中2 組的野生型和突變型純合行各10 行, 每行各10 株完整植株進行考種, 在數據分析時刪除相應性狀的最大值和最小值, 取平均值進行t檢驗, 若P<0.05 則代表此性狀有顯著差異。在北京順義夏季長日照條件下大田種植, 黃化性狀與短日照相同, 因此進行同樣的考種和數據分析, 并對2 種條件下獲得的數據進行對比測驗。

1.3 BSA-Seq 候選區間定位方法

根據2 個親本及F2:3葉片表型, 構建2 個親本和2 個葉色差異表型混池, 親本池分別為10 株野生型(ofc)與10 株突變型(ly); 子代池在F2:3植株中,選擇不同黃葉純合行每行各選1 株, 混合51 單株個體構建突變型混池, 不同綠葉純合行每行各選1 株,混合48 單株個體構建野生型混池。具體操作步驟如下: 首先, 取植株頂端幼葉, 利用CTAB 法分別提取單株葉片DNA; 檢測DNA 質量和濃度, 將不同池的植株DNA 等量混合構建出4 個混池。由北京百邁克生物科技有限公司測序。利用 Illunima Casava 1.8 進行堿基識別分析, 采用雙端150 bp 測序策略進行基因組測序。親本池測序深度為10×, 子代混池測序深度為50×。參考基因組為Wm82.a2.v1版本的大豆基因組[39]。使用GATK[40]軟件工具來檢測SNP, 利用SnpEff[41]軟件進行變異注釋和預測變異影響。

1.3.1 ED 關聯分析 歐式距離(Euclidean Distance, ED)算法, 是利用測序數據尋找混池間存在顯著差異標記, 并以此評估與性狀關聯區域的方法。理論上, 野生型混池和突變型混池之間除了目標性狀相關位點存在差異, 其他位點均趨向于一致, 因此非目標位點的ED 值應趨向于0, ED 值越大表明該標記在兩混池間的差異越大, ED 方法的計算公式如下所示, 式中mut 表示堿基在突變混池中的頻率, wt表示堿基在野生型混池中的頻率。

1.3.2 index 關聯分析 index 是通過混池間的基因型頻率差異進行標記關聯分析的方法, 尋找混池之間基因型頻率的顯著差異并用Δindex 統計。計算公式如下所示, 式中aa 為子代的野生池即綠葉混池,ab 代表子代的突變池即黃葉混池, M 和P 分別為野生型親本的等位基因和突變型親本的等位基因在各自池中出現的read 數目。通過Δindex 可以觀察每個位點在突變型和野生型混池之間的差異。標記與性狀關聯度越強, Δindex 越接近于1。

為檢測SNP 和InDel 的可靠性, 對相應位點進行驗證。首先在Soybase (http://www.soybase.org/)數據庫中查找參考基因組中的變異位點上下游特異序列, 通過DNAMAN 軟件(https://www.lynnon.com/)設計引物, 對PCR 擴增產物進行測序, 鑒定目標位點是否存在。引物設置如表1。

表1 7 個變異位點上下游引物序列Table 1 Upstream and downstream primer sequences of 7 mutation sites

1.4 圖位克隆

檢測SNP/InDel 位點準確性后, 開發分子標記。采用CTAB 法提取對應單株葉片的DNA, 以在F2:3群體中的葉色為表型開展基因型檢測, 根據基因型和表型的關系進行圖位克隆, 再將標記的基因型與表型進行單標記卡方檢測分析, 分析公式為

式中,O為表型與基因型相同的綠葉或黃葉行數;E為鑒定出基因型的綠葉或黃葉行數

1.4.1 競爭性等位基因特異性PCR (kompetitive allelespecific PCR, KASP) 每個標記設計2 條SNP 特異性引物(FAM/HEX)和一條通用引物Rev,FAM 尾部添加能夠與FAM 熒光結合的特異性序列(5′-GAAGGTGACCAAGTTCATGCT-3′), F-HEX 尾部添加能夠與 H E X 熒光結合的特異性序列(5′-GAAGGTCGGAGTCAACGGATT-3′), 引物序列見表2, KASP 擴增反應體系為10 μL, 具體為: 5 μL 2 × KASP Master Mix, 0.14 μL KASP Primer Mix, 1.0 μL 模板DNA (4～50 ng μL–1), 3.86 μL ddH2O。反應程序為: 95℃熱處理10 min; 95℃變性20 s, 61～55℃退火和延伸60 s, 10 個循環(每循環降低0.6℃); 95℃變性20 s, 55℃退火和延伸40 s, 28 個循環; 4℃避光保存, 經熒光信號采集后即可基因分型。

表2 定位區間內7 個標記的引物序列Table 2 Primer sequences of seven markers located in the interval

1.4.2 衍生的酶切擴增多態性標記(derived cleaved amplified polymorphic sequence, dCAPS) 選擇在InDel 位點(AGA/A)處設計 dCAPS 標記, 引物為241500D-F 和241500D-MobI-R2, 命名為241500D。下游引物241500D-MboI-R2 序列中引入1 個錯配堿基G, 使其在該位點形成了MboI (GATC)的酶切位點。在綠色基因型(caaAGAta)材料中, PCR 擴增產物可被MboI 酶切成155 bp 和28 bp 片段, 黃色基因型(caaAta)材料中, 因不含有MboI 酶切位點而不能被酶切。對PCR 擴增產物酶切, 經聚丙烯酰胺凝膠電泳即可基因分型。

1.5 RNA-Seq 差異基因分析方法

根據2 個親本及F3:4葉片表型建池, 在鼓粒后期取頂部新鮮葉片, 親本池分別為10 株野生型(ofc)與10 株突變型(ly); 在F3:4群體中, 子代池選擇純合行黃葉10 個單株構建突變型混池, 綠葉10 個單株構建野生型混池。利用TRIzol Kit 提取葉片總RNA, 采用Truseq RNA 試劑盒進行mRNA 純化并構建cDNA文庫, 對獲得的cDNA 文庫進行PCR 擴增富集, 采用分子生物學先進設備檢測RNA 樣品的純度、濃度和完整性, 以保障使用合格的樣品進行轉錄組測序。庫檢合格后, 不同文庫按照目標下機數據量進行pooling, 用Illumina 平臺進行測序。RNA 提取和測序由百邁克生物科技有限公司完成。經過質量控制之后得到的高質量的Clean Data。再利用Tophat將Clean reads 比對到大豆參考基因組, 允許有2 個堿基的錯配[42-43]。

1.5.1 變異位點分析 基于各樣品 reads 與參考基因組序列的Hisat2 比對結果, 使用GATK[45]軟件識別測序樣品與參考基因組間的單堿基錯配, 識別潛在的 SNP 或InDel 位點。在不同葉色池間篩選變異位點, 再在相同葉色池間確定相同的變異位點,最終獲得可靠的變異位點。

1.5.2 差異表達分析 差異表達分析得到的基因集合叫做差異表達基因集, 使用“A_vs_B”的方式命名。根據兩組樣品之間表達水平的相對高低,差異表達基因劃分為上調基因(up-regulated gene)和下調基因(down-regulated gene)。上調基因在樣品(組) A 中的表達水平高于樣品(組) B 中的表達水平;反之為下調基因, 無表達差異基因為(no-differential genes)。表達定量的結果以FPKM(fragments per kilobase ofexon model per million mapped fragments) 為單位, 在得到差異檢驗的FDR 值同時, 根據基因的表達量(FPKM 值)計算該基因在不同樣本間的差異表達倍數。以“FDR≤0.001 和|log2Ratio (Fold Change)|≥1”作為閾值來判斷基因差異表達的顯著性[44]。根據不同池得到的差異表達基因, 在野生型親本vs 突變型親本、野生型子代vs 突變型子代、野生型親本vs 突變型子代、突變型親本 vs 野生型子代中篩選差異表達基因,這4 個結果取交集得到不同葉色間差異基因; 在野生型親本vs 野生型子代、突變型親本vs 突變型子中的差異表達基因, 這2 個結果都不是不同葉色間的差異表達基因, 因此應取并集, 并在取交集得到的基因中排除, 即得到野生型和突變型間準確的差異表達基因。

2 結果與分析

2.1 大豆鼓粒后期葉片黃化性狀分析

2.1.1 遺傳分析 鳳交66-12、冀豆12、鄭1307和ofc分別與ly雜交, F1葉色表現為綠色, 表明此性狀為隱性性狀。ofc×ly群體F2:3的大豆鼓粒后期葉色的綠色純合行、黃色純合行和雜合行行數分別為212、203 和418, 對其進行適合性測驗, χ2(1:2:1)=0.213< χ20.05,2=5.99,P=0.9, 說明性狀符合1∶2∶1 的概率為90%, 黃化性狀由1 對隱性基因控制。

2.1.2 黃化性狀的變化特點 觀察F2:3鼓粒后期不同植株、葉、莢和種子的黃化性狀, 從表型開始出現差異到葉片凋落, 顏色有明顯差異; 其中突變型為葉脈間開始失綠變黃到整個植株失綠變黃, 而后葉片褐變凋落, 野生型為正常變黃凋落(圖2-A),此表型是一個漸變的過程, 突變型出現黃色表型早;從衰老情況看, 突變型葉片更容易褐化凋落(圖2-B);從突變型的成熟莢和種子上看, 突變型比野生型發白(圖2-C, D)。

圖2 野生型與突變型鼓粒后期植株整體的差異變化(A)、葉色變化差異(B)、莢色變化差異(C)和成熟種子的差異(D)Fig.2 Differences in overall changes of plants in soybean filling later period of wild type and mutant tympanum (A), leaf color change (B), pod color change (C), and mature seeds (D)

大豆鼓粒后期水分和干物質的測定結果表示(表3), 突變型與野生型水分差異不明顯, 但突變型的干重較高, 結合植株成熟后的考種數據, 認為此性狀使鼓粒時的干物質提前積累。

表3 突變型和野生型鼓粒后期水分和干物質的差異Table 3 Differences in moisture and dry matter at the late stages of mutant and wild tympanums

對大豆鼓粒后期的葉片SPAD 值進行分析(圖3)后, 發現至表型出現后隨著時間的推移, 二者SPAD值差異越發明顯, 直到兩表型均測不出SPAD 值為止, 二者的SPAD 值均有明顯差異, 說明葉綠素含量在二者間存在顯著差異。

圖3 大豆鼓粒后期葉片SPAD 值的變化Fig.3 Changes of SPAD values of leaves at the late stage of soybean drum grain

2.1.3 其他農藝性狀比較 對矮皺葉和高展葉兩組的野生型和突變型的考種數據分析(表4)認為短日照條件下差異性狀較多, 矮皺葉組突變型株高、莢數、粒數和百粒重都顯著低于野生型, 高展葉組突變型株高、莢數、粒數和百粒重也顯著低于野生型, 分枝數顯著高于野生型; 長日照下差異性狀較少, 矮皺葉組只有莢數突變型顯著高于野生型, 高展葉組突變型分支和粒數顯著高于野生型, 百粒重顯著低于野生型。

表4 野生型和突變型的考種數據分析Table 4 Analysis of wild-type and mutant-type agronomic trait data

2.2 基于BSA-Seq 篩選黃化性狀的定位區間

2.2.1 數據質控 2 個親本池與2 個子代池通過Illumina HiSeq 重測序, 測序共得到的Clean Data 為141.04 Gbp, Q30 達到90.56。樣品與參考基因組平均比對效率為98.70%, 平均覆蓋深度為30.50X, 基因組覆蓋度為99.33% (至少一個堿基覆蓋), 樣品AT、CG 堿基基本不發生分離且各插入片段大小的分布呈單峰的正態分布, 說明測序數據質量合格, 與大豆參考基因組比對效率較高, 可用于后續分析。

2.2.2 關聯分析 進一步過濾SNP/InDel, 得到高質量可信SNP/InDel 位點12,160/3218 個。采用ED 算法和index 2 種算法對這些高質量SNP/InDel進行關聯。利用兩混池間基因型存在差異的SNP/InDel 位點, 統計各個堿基在不同混池中的深度, 并計算每個位點ED 值, 為消除背景噪音, 對原始ED 值進行5 次方處理[45], 然后采用DISTANCE方法對 ED 值進行擬合, 取所有位點擬合值的median+3SD 作為分析的關聯閾值[46], 計算得0.50/2.09。根據關聯閾值判定, SNP/InDel 的ED 關聯結果如圖4-A 所示, index 關聯結果如圖4-B 所示。

圖4 SNP/InDel 的ED 和index 關聯分析結果Fig.4 ED and index correlation analysis of SNP/InDel

根據SNP 的關聯分析在19 號染色體定位區間4.96 Mb, 包括636 個基因。根據InDel 關聯分析定位了包括19 號染色體在內的8 個區間, 2 種關聯分析方法得到的結果取交集, 僅有1 個區域(表5), 總長度2.23 Mb, 共包含277 個基因。將測序分析結果和關聯分析結果可視化, 并將樣品的變異結果及BSA 關聯分析結果使用circus 軟件(http://circos.ca/)作圖(圖5)。

表5 SNP/InDel 分析結果利用2 種方法在19 號染色體獲得關聯區域Table 5 Associated region on chromosome 19 based on SNP/InDel analysis by two methods

圖5 SNP (A)和InDel (B)在染色體上結果的可視化分布Fig.5 Visual distribution of SNP (A) and InDel (B) results on chromosomes

2.3 開發標記加密定位區間

根據BSA-Seq 分析得到的突變位點, 利用已成功開發的6 個Kasp 標記和一個dCaps 標記對F2:3純合行行混進行標記-表型連鎖分析(圖6), 將區間從2.23 Mb 縮短為物理位置為47,192,234～48,945,702,一共1.75 Mb, 其中含有219 個基因, 其中有SNP 突變的基因5 個。結合單標記分析方法(表6)結果顯示,Glyma.19G223400到Glyma.19G237000之間標記均是χ2< χ20.05,1=3.84, 其中標記237,000 的P值最大,說明在223,400 到237,000 之間的標記基因型與表型相符概率大, 認為其與目標基因的連鎖程度最高。

表6 7 個標記在F2:3 純合行的單標記分析Table 6 Single label analysis of seven markers in F2:3 homozygous rows

圖6 7 個標記在F2:3 純合行的表型連鎖分析Fig.6 Phenotypic linkage analysis of seven labeled homozygous rows at F2:3

2.4 RNA-Seq 差異基因

2.4.1 數據質控 經過測序質量控制, 共得到24.81 Gb Clean Data, 各樣品Q30 堿基百分比均不小于 94.94%。樣品與參考基因組從比對結果統計來看,各樣品的 Reads 與參考基因組的比對效率在52.50%～96.53%之間。通過檢驗插入片段在基因上的分布, reads 在所有表達的基因上的分布呈現均一化分布, 在插入片段大小的分布呈單峰的正態分布,使用各樣品的Mapped Data 對檢測到的不同表達情況的基因數目飽和情況進行模擬, 隨著測序數據量的增加, 數值越趨近于1, 這說明mRNA 片段化的隨機性高、無mRNA 降解情況, 插入片段長度的離散程度高, 文庫容量和Mapped Data 充足[47]?？捎糜诤罄m候選區間的基因篩選。

2.4.2 變異位點分析 對RNA-Seq 得到的變異基因進行篩選, 在已得區間內有4 個SNP 差異的基因(表7 和表8),Glyma.19G223400、Glyma.19G233000和Glyma.19G237000均發生在UTR 區域, 可能導致RNA 二級結構發生改變, 進而影響RNA 的穩定性和翻譯效率,Glyma.19G236800則是在CDS 區發生改變, 使第185 位脯氨酸變成亮氨酸, 改變了蛋白的鋅指結構域, 從而影響蛋白的功能。

表7 SNP 突變的4 個候選基因Table 7 Four candidate genes for SNP mutations

表8 4 個SNP 基因的表達信息Table 8 Relative expression information of four SNP genes

2.4.3 差異表達基因分析 通過不同葉色池間差異表達基因取交集, 得到2948 個基因(圖7-A); 相同葉色池間差異表達基因取并集, 得到4780 個基因(圖7-B), 而后從取交集的基因中排除取并集中的相同基因, 得到2392 個與大豆鼓粒后期葉色差異相關基因, 這些差異表達基因與獲得的區間取交集, 得到8 個基因(表9)。

表9 區間內8 個差異表達基因Table 9 Eight differentially expressed genes in the interval

圖7 野生型和突變型的親本或子代間的差異基因取交集(A)親本和子代的野生型或突變型的差異基因取并集(B)Fig.7 Intersection of differential genes between parents or offspring of wild-type and mutant type (A) combination of differentiating genes of wild-type or mutant type of parents and offspring (B)

2.5 候選基因生物信息學分析

2.5.1 表達譜和蛋白序列分析 為了進一步探究共定位區間內候選基因在各個組織中的表達情況,利用phytozome (/phytozome.jgi.doe.gov/)數據庫查詢了4 個SNP 差異基因和8 個表達差異基因在Williams 82中的表達數據, 從構建不同部位表達圖譜(圖8)結合蛋白功能注釋(表10)中可以看出,Glyma.19G237000、Glyma.19G225900、Glyma.19G226800和Glyma.19G228200在各部位中的表達量都低, 它們與核苷三磷酸水解和碳水化合物代謝過程相關;Glyma.19G222400、Glyma.19G223400和Glyma.19G233000表達模式相近, 均是花中表達低其他部位表達高, 其與乙酰轉移和信號轉導有關,Glyma.19G224500與這些基因的表達模式相反, 在花中表達低其他部位高, 同樣也與乙酰轉移可能有關;Glyma.19G230600在芽頂端分生組織、莖和種子中表達相對較高, 其與蛋白磷酸化有關;Glyma.19G228300在莖和果莢中表達較高,Glyma.19G236800在葉以及根、根瘤和根毛中表達較高, 與泛素化蛋白降解有關;Glyma.19G241200在花、葉、果莢以及芽頂端分生組織中表達較高, 與?；D移相關。

表10 候選基因蛋白序列BLAST 分析Table 10 BLAST analysis of candidate gene protein sequences

圖8 12 個候選基因的表達圖譜Fig.8 Relative expression profile of 12 candidate genes

3 討論

3.1 水分和干物質的測定方法

有試驗表明, 利用莢果厚度可以模擬大豆的鼓粒進程[48], 大豆莢果發育過程有2 個部分重疊的階段: 結莢(建庫)階段和鼓粒(籽粒充實)階段。莢長和莢寬的增加主要集中在結莢階段, 莢果長度和寬度達到最大值時大豆的干重只有最終干重的4%[49]。莢果厚度增加略晚于莢長和莢寬, 而與大豆鼓粒相一致[50]。這說明莢果可以在一定程度上反應鼓粒進程, 而本試驗的性狀在鼓粒后期, 也就是粒鼓滿后表型開始發生變化, 因此在試驗方法上, 選擇使用莢果的鮮重來模擬進程, 可以消除生長上的差異,從而探究水分和干物質的變化。

3.2 大豆鼓粒期干重變化

鼓粒期之前, 豆莢發育從葉片獲得較多的光合產物來構建自身, 并通過較強的呼吸作用為其自身發育提供能量。而進入鼓粒期后, 葉片光合能力達到最大[51-53], 在鼓粒前期, 豆莢具有較強的呼吸作用, 呼吸底物除豆莢自身光合形成的光合產物外,主要來源于葉片光合產物, 此時期豆莢的光合產物對籽粒重的貢獻較大(貢獻率為7.34%～15.06%)。在鼓粒中期(開花26～40 d), 葉片干重先上升后再下降,光合產物主要向籽粒輸送, 此時期籽粒干重增加迅速, 鼓粒速度快, 粒重增加也快, 以生殖生長為主;在鼓粒末期(開花40 d 后), 葉片干重下降, 鼓粒速度下降, 但粒重仍在增加, 單株籽粒產量呈上升趨勢,營養體生長停止, 生殖生長還在繼續[54]; 本實驗研究的性狀發生在鼓粒后期R5 至R6 期, 也就是在籽粒生殖生長階段, 此時的SPAD 值說明黃化性狀的葉綠素含量下降, 認為在營養體生長停止后野生型的植株依然繼續一段生殖生長, 而突變型的植株不再進行生殖生長直接將干物質運往籽粒, 因此導致百粒重上發生顯著的差異。雖然有報道認為百粒重與籽粒含水量關系密切, 呈顯著負相關[54], 但在本研究中, 鼓粒后期黃化與正常性狀的水分差異不明顯, 也就說明了黃化性狀的葉片在營養體生長停止后, 干物質迅速運往籽粒, 但生殖生長階段的葉片光合時間變短, 最終運往籽粒的總干物質減少, 導致百粒重下降。因此本研究的黃化性狀是犧牲籽粒部分干重為代價, 提前將干物質積累完畢?？梢詰糜诠牧：笃趹獙夂蜃兓任粗蛩貢r促使葉片干物質提前轉移完畢, 相對提高一些產量。

3.2 候選基因預測

大豆鼓粒期可溶性糖含量的變化是植物體內碳水化合物代謝的重要指標, 既能夠反映碳水化合物的合成情況, 也能夠說明碳水化合物在植物體內的運輸情況[55-57]。大豆鼓粒期葉片和莢皮都是糖代謝相對比較旺盛的部位, 為籽粒運輸更多的碳水化合物奠定基礎, 可溶性糖含量為葉片>莢皮>莖稈>根系[58]。本試驗對RNA-Seq 得到的2392 個葉色差異表達基因進行生物信息學COG 和KOG 的功能分類分析, 發現了4 個碳水化合物的運輸和代謝功能相關基因, 為Glyma.19G233000、Glyma.19G222400、Glyma.19G224500和Glyma.19G225900, 它們有可能是大豆鼓粒期后期黃化性狀的目的基因; 在KEGG通路富集和KEGG 分類分析中認為與植物激素信號轉導相關的候選基因Glyma.19G223400、Glyma.19G236800和Glyma.19G230600也有可能是大豆鼓粒期后期黃化性狀的目的基因, GO 富集分析認為跟質膜組分相關, 而Glyma.19G236800也同樣是質膜組分, 且4 個SNP 基因分析發現, 3 個UTR 區域突變并沒有改變表達量,Glyma.19G236800有著蛋白序列的變異, 因此認為其更可能是黃化性狀的目的基因。

4 結論

突變體ly的鼓粒后期黃化性狀受單個隱性核基因控制, 基于BSA-Seq、RNA-Seq 和圖位克隆將該基因定位在19 號染色體上, 獲得1.84 Mb 的區間,區間內有219 個基因, 有SNP 的基因4 個, 有差異表達基因8 個, 它們可能與大豆鼓粒后期葉片黃化衰老性狀相關, 但還需后期試驗篩選和驗證這12 個候選基因。本研究為探索大豆葉片黃化的分子調控機制奠定了基礎。