大豆葉型性狀全基因組關聯分析與候選基因鑒定

王 瓊 朱宇翔,2 周密密 張 威 張紅梅 陳 新 陳華濤,*崔曉艷,*

1 江蘇省農業科學院經濟作物研究所 / 江蘇省高效園藝作物遺傳改良重點實驗室, 江蘇南京 210014; 2 揚州大學園藝園林學院, 江蘇揚州 225009

大豆是世界重要的經濟作物, 是人類和家養動物蛋白質和油脂的重要來源[1]。近年來, 我國大豆對外依存度較高, 供需形勢日益嚴峻[2]。單產水平低是我國大豆生產面臨的最大困境, 培育高產大豆品種、提高大豆單產成為扭轉我國大豆被動局面的首要任務[3]。作物的產量是復雜性狀的集合, 受許多形態、生理和農藝性狀的影響。葉片是植物進行光合作用的主要器官, 合適的葉片形態有利于塑造理想株型, 增強光合能力, 從而提高作物的產量[4]。

理想株型育種是實現大豆超高產的重要途徑[5]。葉型是大豆株型性狀的重要組成部分, 影響冠層結構、光合作用效率, 并最終影響產量[6]。通過對葉型結構進行優化, 可以實現培育高產大豆新品種的目標。彭玉華等[7]認為, 按葉型垂直分布不同, 大豆的株型可分為3 種類型, 即下闊葉-上闊葉、下闊葉-上窄葉和下窄葉-上窄葉, 其中以下窄葉-上窄葉類型的單莢粒數最多, 以下闊葉-上闊葉和下闊葉-上窄葉類型的百粒重最大。盡管下闊葉-上窄葉類型的單莢粒數略低于下窄葉-上窄葉類型, 但卻顯著高于下闊葉-上闊葉類型, 因此, 以下闊葉-上窄葉為目標進行選擇, 可以在穩定其他產量因素的前提下, 使育成品種的單莢粒數顯著提高, 從而實現對大豆品種的產量改良。

然而, 目前大豆葉型相關的研究大多集中在單一類型的葉形性狀上, 關于葉型垂直分布的異形葉的研究還較少。Chen 等[8]利用由葉片形狀明顯不同的3 個大豆材料衍生構建的2 個重組自交系(recombinant inbred lines, RIL), 通過連鎖分析共檢測到8 個與異形葉性狀相關的QTL, 再根據2 個作圖群體3 個親本的序列差異、RT-qPCR 分析和基因功能注釋分析,鑒定出4 個異形葉候選基因, 其中包含大豆葉形調控基因Ln。

近年來, 隨著高通量測序和高密度SNP 基因分型等技術的快速發展, 全基因組關聯分析(Genomewide association study, GWAS)已廣泛應用于大豆復雜性狀的遺傳解析[9-12]。本研究利用283 份大豆種質資源數據, 通過分析一系列葉形相關性狀, 包括葉長(leaf length, LL)、葉寬(leaf width, LW)、葉形指數(leaf shape, LS)和異形葉指數(heterophylly index, HI)等, 明確異形葉對大豆產量形成的影響; 檢測與目標性狀關聯的區域與位點, 探究大豆葉型相關性狀形成的分子機制, 包括葉片的大小、形狀及其層次分布等。本研究將有助于解析大豆理想株型的組成要素, 挖掘鑒定與大豆葉型性狀關聯的等位變異和候選基因, 以期為大豆產量改良提供指導。

1 材料與方法

1.1 試驗材料與表型鑒定

供試群體為283 份大豆種質資源組成的自然群體, 包含野生種質19 份、地方品種52 份和212 份栽培品種, 這些大豆材料代表了豐富的遺傳多樣性。所有材料于2021 年6 月25 日和2022 年6 月20 日在江蘇省農業科學院六合基地種植2 年。種植3 個重復, 隨機區組設計, 大田常規管理。

測量上部和下部葉片的葉長、葉寬, 其中上部葉測量倒三葉, 下部葉測量第一復葉, 分別測量3 個單葉(圖1)。用上部葉長除以上部葉寬得出上部葉形指數, 用下部葉長除以下部葉寬得到下部葉形指數, 再計算上下部葉形指數的比值, 得到異形葉指數。

圖1 試驗設計示意圖Fig. 1 Schematic representation of the design in the experiment

1.2 重測序與基因分型

利用CTAB 法[13]提取全部供試群體材料的葉片基因組DNA, 利用Illumina HiSeq 高通量測序儀進行測序, 每個樣品的平均測序數據量不低于10 倍基因組覆蓋深度。使用Trimmomatic 軟件[14]對測序產生的原始數據進行過濾, 去除接頭污染和低質量數據后, 將得到的高質量數據比對大豆參考基因組,檢測供試群體SNP 標記。具體SNP 檢測流程參考本課題組前期相關工作[15]。

1.3 全基因組關聯分析

篩選最小等位基因型頻率(minor allele frequency,MAF)≥5%、缺失率＜10%的SNP 位點, 并利用Beagle軟件[16]進行基因型填充后, 結合7 個葉型相關性狀的表型數據, 運用EMMAX 軟件[17]中的混合線性模型(linear mixed model, MLM)進行GWAS 分析。

1.4 候選基因鑒定

根據大豆群體連鎖不平衡衰減距離, 選取顯著SNP 位點及其上下游一定區段(120 kb)作為該位點候選基因分析區間, 結合基因在大豆各個組織的表達量數據、基因的注釋信息和已報道相關QTL 的連鎖定位, 挖掘鑒定與葉型性狀關聯的遺傳標記或候選基因。

1.5 實時熒光定量PCR

使用北京天根生化科技有限公司RNAsimple 總RNA 提取試劑盒提取根、根瘤、莖、葉片、莖端分生組織、花、莢和種子的總RNA, 使用PrimeScript RT reagent kit (TaKaRa)合成cDNA。利用Primer Premier 5.0 軟件設計引物, Forward: 5'-GCCCCAGG CTGAAGAAAATC-3', Reverse: 5'-CCACTCGAGGG AAGAAGTAGC-3', 使用Tubulin基因作為內參, 正向引物序列為: 5'-GGAGTTCACAGAGGCAGAG-3',反向引物序列為: 5'-CACTTACGCATCACATAGCA-3'。預混液使用SYBR Green Realtime PCR Master Mix (Toyobo), 采用ABI 7500 Real-Time PCR System分析數據。用2-ΔΔCt法計算基因的相對表達量, 每個樣品3 次重復。

2 結果與分析

2.1 葉型相關性狀表型分析

283 份材料于2021 年、2022 年在江蘇南京連續種植2 年, 分別測量R6～R8 期群體材料上下部葉片的葉長、葉寬, 計算葉形和異形葉指數。對2 個環境的上下部葉長、葉寬、葉形以及異形葉指數等7個葉型相關性狀的表型值進行統計分析發現, 大豆種質資源葉型性狀的變異幅度較大, 單個環境的上部葉長極差范圍為10.57～11.00, 上部葉寬極差范圍為7.20～7.63, 上部葉形極差范圍是1.53～1.89; 下部葉長的極差范圍為8.10～10.30, 下部葉寬極差范圍為5.23～6.70, 下部葉形極差范圍是1.02～1.71, 異形葉指數的極差范圍為1.08～1.21 (表1 和圖2)。7 個性狀的表型值都呈現正態分布, 表明葉型相關性狀為多基因控制的數量性狀(圖2)。2 個環境的性狀間存在顯著相關性(P≤0.001), 表明本研究考察的表型數據具有較高的可靠性和重復性(表2)。

表1 大豆群體葉型相關性狀的統計分析Table 1 Statistical analysis of leaf characteristics traits of soybean association panel

表2 葉型相關性狀的相關系數Table 2 Correlation coefficients of leaf characteristics traits

圖2 大豆關聯群體葉型相關性狀的分布Fig. 2 Distribution of leaf characteristics traits of soybean association panel

2.2 不同葉型性狀相關性分析

對不同葉型性狀進行相關性分析(圖3)發現, 上下部葉片的葉長和葉寬之間都呈顯著正相關, 而葉寬與葉形性狀存在極顯著負相關; 下部葉長與上部葉長之間、下部葉寬與上部葉寬之間均呈顯著正相關, 同時上下部葉形性狀間也呈顯著正相關, 說明大豆種質資源葉形垂直分布以下闊葉—上闊葉和下窄葉—上窄葉類型為主; 此外, 下部葉寬與上部葉形性狀之間、上部葉寬與下部葉形性狀之間均呈顯著負相關, 同樣說明大豆葉形分布以上下都為闊葉或窄葉為主; 值得注意的是, 下部葉長與下部葉形性狀之間, 存在極顯著正相關, 而上部葉長與上部葉形性狀無顯著相關性; 異形葉指數與上部葉形性狀存在顯著正相關, 而與上部葉寬之間、與下部葉長和下部葉形性狀之間均呈顯著負相關。

圖3 大豆關聯群體葉型性狀間的相關系數Fig. 3 Correlation coefficients among leaf characteristics traits in soybean

2.3 全基因組關聯分析

利用3,319,306 個高質量SNP (MAF≥0.05)進行基因分型。結合葉型性狀表型數據, 使用EMMAX軟件的混合模型進行全基因組關聯分析, 合并屬于同一關聯信號的SNP, 共檢測到181 個顯著相關SNP 位點, 其中能夠在2 個環境或多個性狀中重復檢測到的GWAS 信號有18 個, 分布在5、10、13、16、17、19 和20 號染色體上(表3 和附圖1)。位于19 號染色體的上部葉寬關聯信號19_23.9 在2021 和2022 兩個環境中被重復檢測到; 位于5 號染色體的關聯信號在2021 年下部葉長和下部葉寬2 個性狀中被重復檢測到; 位于10 號染色體的關聯信號能夠在2021 年的上部葉形和異形葉指數性狀中被重復檢測;位于13 號染色體上的GWAS 信號在2022 年上部葉寬和2021 年上部葉形性狀中被重復檢測到; 16 號染色體上的關聯信號在2021 年的上下部葉長性狀中被重復檢測; 位于17 號染色體的關聯信號能夠在2022 年上下部葉長性狀中被重復檢測到; 另外, 位于20 號染色體上的GWAS 信號在2022 年上下部的葉形性狀中被重復檢測; 值得注意的是, 位于19 號染色體上的關聯信號19_45.1 在2021 年上部葉形和異形葉指數、2022 年上部葉寬和上部葉形等4 個環境或性狀中被重復檢測到。表明大豆葉片的發育過程, 受到復雜的遺傳網絡調控。

表3 葉型相關性狀的關聯信號Table 3 GWAS signals of leaf characteristics traits

在所有葉形性狀相關的SNP 位點中, 最為顯著的位點為位于2 號染色體上的上部葉寬相關位點(P=6.15E-9)。葉長性狀最顯著相關的 SNP 位點(Chr20:37941224)位于20 號染色體(P=1.49E-7)。對于異形葉指數性狀, 最為顯著相關的SNP 位點位于16 號(P=2.58E-8)和19 號(P=2.84E-8)染色體, 其中位于19 號染色體上的相關SNP (Chr19:45155943)同樣也是葉形指數性狀(2021 上部葉形指數性狀)最顯著相關的位點。綜上所述, 這些顯著位點與特定葉型性狀關聯, 特別是在不同環境中多次檢測到的相關位點, 將有助于闡明大豆中葉型調控和異形葉形成的機制, 為通過葉型設計培育高產大豆品種提供優異的基因資源。

2.4 關鍵候選基因鑒定

利用檢測到的與葉型性狀相關的位點, 通過基因的表達譜數據、擬南芥中同源基因的功能注釋,確定與大豆葉片發育和異形葉形成相關的候選基因?；诖蠖够蚪M注釋信息, 在19 號染色體Chr19:45151266 位點下游11.7 kb, 鑒定到上部葉寬、上部葉形和異形葉指數性狀的候選基因Glyma.19G 194100?；蚬δ茏⑨尡砻?Glyma.19G194100是LITTLE ZIPPER 3 (ZPR3)基因的同源基因, 并且該基因在幼葉中具有較高的表達量。在擬南芥中,ZPR3是莖端分生組織(Shoot apical meristem, SAM)發育過程中HD-ZIPIII 轉錄因子的競爭性抑制基因,該基因通過形成非功能性異二聚體起作用, 參與莖端分生組織的發育和側生器官發生, 對莖端分生組織中干細胞發揮正常功能至關重要[18]。葉片是莖端分生組織產生的第一類側生器官, 在植物的發育中具有重要地位。因此Glyma.19G194100可能是調控大豆葉片形狀和垂直分布的關鍵候選基因。另外, 在19號染色體的相關位點Chr19:45155943 上游109.9 kb發現基因Glyma.19G192700, 該基因是 Growthregulating factor 4 (GRF4)基因的同源基因。生長調節因子GRF 參與調節植物生長發育的各個方面, 包括葉片大小、葉片衰老、應激反應、雌蕊發育和種子大小等[19-24]。GRF4 作為轉錄調節因子GRF 家族的成員, 能夠調節植物的株型和種子大小、氮肥利用效率等[25-26]。擬南芥葉片在近端和遠端軸上具有特定的生長極性, miR396 在幼嫩發育葉片的遠端優先積累并靶向GRF 從而調控這種特定的葉片生長發育模式[27]。利用熒光定量 PCR 分析Glyma.19G 192700在大豆各組織中的表達模式(圖4)發現, 該基因在莖端分生組織中高表達。SNP 注釋分析表明,Glyma.19G192700第3 個外顯子的第168 個堿基發生C→A 替換, 導致編碼的氨基酸從亮氨酸變為異亮氨酸。說明,Glyma.19G192700可能是大豆控制葉型等性狀的關鍵候選基因。

圖4 候選基因Glyma.19G192700 在大豆各組織中的表達模式Fig. 4 Relative expression pattern of candidate genes Glyma.19G192700

3 討論

葉型包括葉片的形狀和垂直分布不僅能夠影響植物對光能的吸收和利用、光合效率, 也可以改變群體的冠層結構和冠層空間光分布, 從而影響通風透氣性和光能利用效率[28-29], 并最終影響產量。本研究通過對283 份大豆種質資源的7 個葉型相關性狀進行連續2 年的考察, 并通過全基因組關聯分析檢測到181 個相關位點。其中能夠在2 個環境或多個性狀中重復檢測到的位點有18 個, 能被重復檢測到的位點較少, 說明大豆葉型性狀可能受環境影響較大。

Tischner 等[30]通過對3 個大豆品種相互雜交得到的3 個RIL 群體進行QTL 分析表明, 大豆的種子敗育、葉形和繁殖等相關性狀由多個顯著數量性狀位點QTL 控制, 占遺傳變異的50%左右。這些QTL分別與影響雌雄不育的Ms1、Ms6和St5基因以及影響小葉數目和葉片形狀的Lf1和ln基因相關。其中,以窄葉命名的ln基因具有多效性, 長期以來ln一直被認為還與大豆的單莢粒數性狀緊密相關[31]。窄葉由隱性基因ln控制, 其顯性等位基因Ln控制闊葉的形成。Ln是擬南芥JAGGED(JAG)基因的同源基因,大豆葉片由闊葉(Ln)到窄葉(ln)的轉變與JAG基因中EAR 基序的氨基酸替換有關[32-33]。本研究在20 號染色體Chr20:36152820 位點上游322.7 kb 發現Ln基因(Glyma.20G116200), 相關性狀為2022 年下部葉形指數(P=7.84E-7)。值得注意的是, 經典的開花抑制因子E1也能夠通過直接抑制葉片發育相關基因TCP14和TCP29參與大豆葉片發育的過程[34]。在本研究中,E1基因(Glyma.06G207800)位于6 號染色體上的相關位點Chr06:20059560 (2021 上部葉寬,P=8.37E-6)的下游, 距離為147.5 kb。Chen 等[8]利用具有異形葉性狀的冀豆17 分別與卵圓形葉的冀豆12和披針形葉的綏農14 構建的2 個RIL 群體, 鑒定到異形葉性狀相關的候選基因4 個, 分別位于1、18、19 和20 號染色體上。其中19 號染色體上的候選基因Glyma.19G194300位于本研究中鑒定的相關位點Chr19:45151266 附近。生長調節因子GRF 具有保守的QLQ 和WRC 結構域, 并且它們中的大多數可以與GRF 互作因子(GRF-Interacting-Factor, GIF)相互作用。此外, 大多數GRF 受miR396 的調控, miR396、GRF 及其相互作用因子GIF 共同組成復雜的調控網絡。在水稻中, OsmiR396c-OsGRF4-OsGIF1 分子調控模塊調節水稻的籽粒大小和產量[25]。本研究中預測的候選基因Glyma.19G192700, 作為擬南芥GRF4基因的同源基因, 可能具有一因多效的作用, 除與異形葉等葉型性狀相關外, 可能也與粒重等性狀相關。

4 結論

本研究初步分析了大豆異形葉等葉型性狀的遺傳特征, 利用全基因組關聯分析檢測到181 個葉型性狀相關SNP 位點, 其中能夠在2 個環境或多個性狀中重復檢測到的位點 18 個。最顯著相關位點Chr19:45155943 位于19 號染色體, 在其附近鑒定到候選基因Glyma.19G192700、Glyma.19G194100, 分別被注釋為 Growth-regulating factor 4 (GRF4)和LITTLE ZIPPER 3 (ZPR3)基因的同源基因。此外, 在其他多個位點附近發現大豆已知葉形和葉片發育相關基因Ln和E1。